
#### ALL ANALYSES BELOW ARE CONSISTENT WITH THE PRE-REGISTRATION ####
#### EXPLORATORY ANALYSES ARE CLEARLY IDENTIFIED AS SUCH 

#### LOAD PACKAGES AND SCRIPTS ####

library(psych) # basic psych package
library(RSA) # plotting of response surfaces
library(multicon) # relationships among multivariate constructs i.e. splithalf reliabilities
library(olsrr) # measure of collinearity for OLS regressions
library(lmtest) # tests assumptions of OLS
library(dplyr)
library(lavaan)

#### LOAD AND PREPARE DATA ####

#### The materials, data, and code for analysis have been deposited at Harvard Dataverse: 
#### https://dataverse.harvard.edu/dataset.xhtml?persistentId=doi:10.7910/DVN/BSUNGB. 
#### To download the csv versions of the data that are required by the code, 
#### researchers need to select the "Original File Format (Comma Separated Values)" option in Dataverse. 
#### The code for analysis will not work if researchers download the tab versions of the data files instead.

#### IMPORTANT!!! set working directory here using setwd command
setwd("~/Dropbox/Emotions in Organizations/Accuracy of Self-Insight/Nature Human Behavior/Submitted Stage 2 Article/Dataverse Documents/")

#### load the data from the initial survey
selfinsight <- read.csv("Stage 2 - selfinsight.csv", header=TRUE, stringsAsFactors = FALSE)
#### note: to analyze the data with the 8th surveys that a few participants complete by mistake, use the "selfinsight - with extra diaries.csv" data file instead

#### calculate average time it took participants to complete the initial survey
meanDuration.initial.survey.seconds <- mean(selfinsight$Duration..in.seconds.)
meanDuration.initial.survey.minutes <- meanDuration.initial.survey.seconds/60
meanDuration.initial.survey.minutes

#### DEMOGRAPHIC VARIABLES ####

# age
selfinsight$age<- as.numeric(selfinsight$age)
selfinsight$age.corrected <- selfinsight$age+12
describe(selfinsight$age.corrected)

# education
table(selfinsight$education)
#1 less than high school
#2 high school graduation or some university
#3 university graduation
#4 master's degree
#5 PhD degree
#6 Professional degree (MD or JD)
#7 Other (please specify)
table(selfinsight$education_7_TEXT)

# income
table(selfinsight$income)
# recode income brackets into values
selfinsight$personalincome <- NA
selfinsight$personalincome[selfinsight$income==1] <- 0
selfinsight$personalincome[selfinsight$income==2] <- 2500
selfinsight$personalincome[selfinsight$income==3] <- 7500
selfinsight$personalincome[selfinsight$income==4] <- 12500
selfinsight$personalincome[selfinsight$income==5] <- 17500
selfinsight$personalincome[selfinsight$income==6] <- 22500
selfinsight$personalincome[selfinsight$income==7] <- 27500
selfinsight$personalincome[selfinsight$income==8] <- 35000
selfinsight$personalincome[selfinsight$income==9] <- 45000
selfinsight$personalincome[selfinsight$income==10] <- 55000
selfinsight$personalincome[selfinsight$income==11] <- 65000
selfinsight$personalincome[selfinsight$income==12] <- 77500
selfinsight$personalincome[selfinsight$income==13] <- 92500
selfinsight$personalincome[selfinsight$income==14] <- 125000
selfinsight$personalincome[selfinsight$income==15] <- 175000
selfinsight$personalincome[selfinsight$income==16] <- 225000
selfinsight$personalincome[selfinsight$income==17] <- 275000
selfinsight$personalincome[selfinsight$income==18] <- 325000
selfinsight$personalincome[selfinsight$income==19] <- 375000
selfinsight$personalincome[selfinsight$income==20] <- 425000
selfinsight$personalincome[selfinsight$income==21] <- 475000
selfinsight$personalincome[selfinsight$income==22] <- 500000
describe(selfinsight$personalincome)
#### create data frames to look at occupations of people with the highest and lowest incomes
selfinsight.lowerincome.temp1 <- subset(selfinsight, income<6); selfinsight.lowerincome <- with(selfinsight.lowerincome.temp1,data.frame(income,occupation),na.rm=TRUE); table(selfinsight.lowerincome$occupation)
selfinsight.higherincome.temp1 <- subset(selfinsight, income>13); selfinsight.higherincome <- with(selfinsight.higherincome.temp1,data.frame(income,occupation),na.rm=TRUE); table(selfinsight.higherincome$occupation)

# gender
table(selfinsight$gender)
# recode open-ended responses into gender categories
selfinsight$gendercoded <- NA;
selfinsight$gendercoded[selfinsight$gender=="male"| selfinsight$gender=="Male"|selfinsight$gender=="MALE"|selfinsight$gender=="m"|selfinsight$gender=="M"|selfinsight$gender=="Male "|selfinsight$gender=="male "|selfinsight$gender=="Male." | selfinsight$gender=="Male/Masc"] <- 1;
selfinsight$gendercoded[selfinsight$gender=="f"| selfinsight$gender=="F"|selfinsight$gender=="FEMAIL"|selfinsight$gender=="female"|selfinsight$gender=="Female"|selfinsight$gender=="FEMALE"|selfinsight$gender=="female "|selfinsight$gender=="Female " | selfinsight$gender=="females"|selfinsight$gender=="woman"|selfinsight$gender=="Woman"] <- 2;
selfinsight$gendercoded[selfinsight$gender=="Male/Genderfluid"| selfinsight$gender=="nonbinary"] <- 3;
selfinsight$gendercoded[selfinsight$gender=="34"| selfinsight$gender=="39"| selfinsight$gender=="49"| selfinsight$gender=="69"| selfinsight$gender=="asexual"] <- 4;
table(selfinsight$gendercoded)
# male (1)
# female (2)
# other (3)
# no valid response (4) (plus one more non-valid response that could not be included in the code)

# ethnicity
table(selfinsight$ethnicity)
# recode number of participants who selected each of the following categories (Caucasian, Asian, African American, Hispanic, and Other)
#1 Caucasian
selfinsight$Caucasian <- NA; selfinsight$Caucasian[selfinsight$ethnicity=="1" | selfinsight$ethnicity=="1,2" | selfinsight$ethnicity=="1,2,3" | selfinsight$ethnicity=="1,2,7" | selfinsight$ethnicity=="1,3" | selfinsight$ethnicity=="1,5" | selfinsight$ethnicity=="1,5,6" | selfinsight$ethnicity=="1,6" | selfinsight$ethnicity=="1,6,7" | selfinsight$ethnicity=="1,7" | selfinsight$ethnicity=="1,8"] <- 1; table(selfinsight$Caucasian)
#2 East Asian 
#3 Southeast Asian 
selfinsight$Asian <- NA; selfinsight$Asian[selfinsight$ethnicity=="1,2" | selfinsight$ethnicity=="1,2,3" | selfinsight$ethnicity=="1,2,7" | selfinsight$ethnicity=="1,3" | selfinsight$ethnicity=="2" | selfinsight$ethnicity=="2,3" | selfinsight$ethnicity=="3" | selfinsight$ethnicity=="3,8"] <- 1; table(selfinsight$Asian)
#6 African American 
selfinsight$AfricanAmerican <- NA; selfinsight$AfricanAmerican[selfinsight$ethnicity=="1,5,6" | selfinsight$ethnicity=="1,6" | selfinsight$ethnicity=="1,6,7" | selfinsight$ethnicity=="4,6" | selfinsight$ethnicity=="4,6,7" | selfinsight$ethnicity=="6" | selfinsight$ethnicity=="6,7" | selfinsight$ethnicity=="6,8"] <- 1; table(selfinsight$AfricanAmerican)
#7 Hispanic
selfinsight$Hispanic <- NA; selfinsight$Hispanic[selfinsight$ethnicity=="1,2,7" | selfinsight$ethnicity=="1,6,7" | selfinsight$ethnicity=="1,7" | selfinsight$ethnicity=="4,6,7" | selfinsight$ethnicity=="4,7" | selfinsight$ethnicity=="6,7" | selfinsight$ethnicity=="7"] <- 1; table(selfinsight$Hispanic)
#4 West Indian 
#5 Middle Eastern 
#8 Other (specify) 
selfinsight$Other <- NA; selfinsight$Other[selfinsight$ethnicity=="1,8" | selfinsight$ethnicity=="3,8" | selfinsight$ethnicity=="6,8" | selfinsight$ethnicity=="8" | selfinsight$ethnicity=="1,5" | selfinsight$ethnicity=="1,5,6" | selfinsight$ethnicity=="4" | selfinsight$ethnicity=="4,6" | selfinsight$ethnicity=="4,6,7" | selfinsight$ethnicity=="4,7" | selfinsight$ethnicity=="5"] <- 1; table(selfinsight$Other)
table(selfinsight$ethnicity_8_TEXT)

# marital status
table(selfinsight$marital)
# Single, never married  (1) 
# In a relationship, but not cohabitating  (2) 
# In a relationship and cohabitating, but not married  (3) 
# Married  (4) 
# Widowed  (5) 
# Divorced  (6) 
# Separated  (7)

#### SCORE TESTS OF ABILITIES ####

#### score the test of emotional abilities
#### response options:
#### Anger  (1) 
#### Disgust  (2) 
#### Fear  (3) 
#### Happiness  (4) 
#### Sadness  (5) 
#### Shame  (6) 
selfinsight$emorecogitem1 <- NA; selfinsight$emorecogitem1[selfinsight$eraitem1==6] <- 1; selfinsight$emorecogitem1[selfinsight$eraitem1!=6] <- 0
selfinsight$emorecogitem2 <- NA; selfinsight$emorecogitem2[selfinsight$eraitem2==3] <- 1; selfinsight$emorecogitem2[selfinsight$eraitem2!=3] <- 0
selfinsight$emorecogitem3 <- NA; selfinsight$emorecogitem3[selfinsight$eraitem3==2] <- 1; selfinsight$emorecogitem3[selfinsight$eraitem3!=2] <- 0
selfinsight$emorecogitem4 <- NA; selfinsight$emorecogitem4[selfinsight$eraitem4==6] <- 1; selfinsight$emorecogitem4[selfinsight$eraitem4!=6] <- 0
selfinsight$emorecogitem5 <- NA; selfinsight$emorecogitem5[selfinsight$eraitem5==1] <- 1; selfinsight$emorecogitem5[selfinsight$eraitem5!=1] <- 0
selfinsight$emorecogitem6 <- NA; selfinsight$emorecogitem6[selfinsight$eraitem6==4] <- 1; selfinsight$emorecogitem6[selfinsight$eraitem6!=4] <- 0
selfinsight$emorecogitem7 <- NA; selfinsight$emorecogitem7[selfinsight$eraitem7==6] <- 1; selfinsight$emorecogitem7[selfinsight$eraitem7!=6] <- 0
selfinsight$emorecogitem8 <- NA; selfinsight$emorecogitem8[selfinsight$eraitem8==1] <- 1; selfinsight$emorecogitem8[selfinsight$eraitem8!=1] <- 0
selfinsight$emorecogitem9 <- NA; selfinsight$emorecogitem9[selfinsight$eraitem9==6] <- 1; selfinsight$emorecogitem9[selfinsight$eraitem9!=6] <- 0
selfinsight$emorecogitem10 <- NA; selfinsight$emorecogitem10[selfinsight$eraitem10==5] <- 1; selfinsight$emorecogitem10[selfinsight$eraitem10!=5] <- 0
selfinsight$emorecogitem11 <- NA; selfinsight$emorecogitem11[selfinsight$eraitem11==4] <- 1; selfinsight$emorecogitem11[selfinsight$eraitem11!=4] <- 0
selfinsight$emorecogitem12 <- NA; selfinsight$emorecogitem12[selfinsight$eraitem12==3] <- 1; selfinsight$emorecogitem12[selfinsight$eraitem12!=3] <- 0
selfinsight$emorecogitem13 <- NA; selfinsight$emorecogitem13[selfinsight$eraitem13==4] <- 1; selfinsight$emorecogitem13[selfinsight$eraitem13!=4] <- 0
selfinsight$emorecogitem14 <- NA; selfinsight$emorecogitem14[selfinsight$eraitem14==1] <- 1; selfinsight$emorecogitem14[selfinsight$eraitem14!=1] <- 0
selfinsight$emorecogitem15 <- NA; selfinsight$emorecogitem15[selfinsight$eraitem15==3] <- 1; selfinsight$emorecogitem15[selfinsight$eraitem15!=3] <- 0
selfinsight$emorecogitem16 <- NA; selfinsight$emorecogitem16[selfinsight$eraitem16==1] <- 1; selfinsight$emorecogitem16[selfinsight$eraitem16!=1] <- 0
selfinsight$emorecogitem17 <- NA; selfinsight$emorecogitem17[selfinsight$eraitem17==2] <- 1; selfinsight$emorecogitem17[selfinsight$eraitem17!=2] <- 0
selfinsight$emorecogitem18 <- NA; selfinsight$emorecogitem18[selfinsight$eraitem18==5] <- 1; selfinsight$emorecogitem18[selfinsight$eraitem18!=5] <- 0
selfinsight$emorecogitem19 <- NA; selfinsight$emorecogitem19[selfinsight$eraitem19==1] <- 1; selfinsight$emorecogitem19[selfinsight$eraitem19!=1] <- 0
selfinsight$emorecogitem20 <- NA; selfinsight$emorecogitem20[selfinsight$eraitem20==5] <- 1; selfinsight$emorecogitem20[selfinsight$eraitem20!=5] <- 0
selfinsight$emorecogitem21 <- NA; selfinsight$emorecogitem21[selfinsight$eraitem21==4] <- 1; selfinsight$emorecogitem21[selfinsight$eraitem21!=4] <- 0
selfinsight$emorecogitem22 <- NA; selfinsight$emorecogitem22[selfinsight$eraitem22==5] <- 1; selfinsight$emorecogitem22[selfinsight$eraitem22!=5] <- 0
selfinsight$emorecogitem23 <- NA; selfinsight$emorecogitem23[selfinsight$eraitem23==3] <- 1; selfinsight$emorecogitem23[selfinsight$eraitem23!=3] <- 0
selfinsight$emorecogitem24 <- NA; selfinsight$emorecogitem24[selfinsight$eraitem24==4] <- 1; selfinsight$emorecogitem24[selfinsight$eraitem24!=4] <- 0
selfinsight$emorecogitem25 <- NA; selfinsight$emorecogitem25[selfinsight$eraitem25==3] <- 1; selfinsight$emorecogitem25[selfinsight$eraitem25!=3] <- 0
selfinsight$emorecogitem26 <- NA; selfinsight$emorecogitem26[selfinsight$eraitem26==5] <- 1; selfinsight$emorecogitem26[selfinsight$eraitem26!=5] <- 0
selfinsight$emorecogitem27 <- NA; selfinsight$emorecogitem27[selfinsight$eraitem27==4] <- 1; selfinsight$emorecogitem27[selfinsight$eraitem27!=4] <- 0
selfinsight$emorecogitem28 <- NA; selfinsight$emorecogitem28[selfinsight$eraitem28==1] <- 1; selfinsight$emorecogitem28[selfinsight$eraitem28!=1] <- 0
selfinsight$emorecogitem29 <- NA; selfinsight$emorecogitem29[selfinsight$eraitem29==6] <- 1; selfinsight$emorecogitem29[selfinsight$eraitem29!=6] <- 0
selfinsight$emorecogitem30 <- NA; selfinsight$emorecogitem30[selfinsight$eraitem30==1] <- 1; selfinsight$emorecogitem30[selfinsight$eraitem30!=1] <- 0
selfinsight$emorecogitem31 <- NA; selfinsight$emorecogitem31[selfinsight$eraitem31==3] <- 1; selfinsight$emorecogitem31[selfinsight$eraitem31!=3] <- 0
selfinsight$emorecogitem32 <- NA; selfinsight$emorecogitem32[selfinsight$eraitem32==2] <- 1; selfinsight$emorecogitem32[selfinsight$eraitem32!=2] <- 0
selfinsight$emorecogitem33 <- NA; selfinsight$emorecogitem33[selfinsight$eraitem33==4] <- 1; selfinsight$emorecogitem33[selfinsight$eraitem33!=4] <- 0
selfinsight$emorecogitem34 <- NA; selfinsight$emorecogitem34[selfinsight$eraitem34==2] <- 1; selfinsight$emorecogitem34[selfinsight$eraitem34!=2] <- 0
selfinsight$emorecogitem35 <- NA; selfinsight$emorecogitem35[selfinsight$eraitem35==5] <- 1; selfinsight$emorecogitem35[selfinsight$eraitem35!=5] <- 0
selfinsight$emorecogitem36 <- NA; selfinsight$emorecogitem36[selfinsight$eraitem36==6] <- 1; selfinsight$emorecogitem36[selfinsight$eraitem36!=6] <- 0
selfinsight$emorecogitem37 <- NA; selfinsight$emorecogitem37[selfinsight$eraitem37==2] <- 1; selfinsight$emorecogitem37[selfinsight$eraitem37!=2] <- 0
selfinsight$emorecogitem38 <- NA; selfinsight$emorecogitem38[selfinsight$eraitem38==4] <- 1; selfinsight$emorecogitem38[selfinsight$eraitem38!=4] <- 0
selfinsight$emorecogitem39 <- NA; selfinsight$emorecogitem39[selfinsight$eraitem39==4] <- 1; selfinsight$emorecogitem39[selfinsight$eraitem39!=4] <- 0
selfinsight$emorecogitem40 <- NA; selfinsight$emorecogitem40[selfinsight$eraitem40==5] <- 1; selfinsight$emorecogitem40[selfinsight$eraitem40!=5] <- 0
selfinsight$emorecogitem41 <- NA; selfinsight$emorecogitem41[selfinsight$eraitem41==2] <- 1; selfinsight$emorecogitem41[selfinsight$eraitem41!=2] <- 0
selfinsight$emorecogitem42 <- NA; selfinsight$emorecogitem42[selfinsight$eraitem42==6] <- 1; selfinsight$emorecogitem42[selfinsight$eraitem42!=6] <- 0
selfinsight$emorecogitem43 <- NA; selfinsight$emorecogitem43[selfinsight$eraitem43==3] <- 1; selfinsight$emorecogitem43[selfinsight$eraitem43!=3] <- 0
selfinsight$emorecogitem44 <- NA; selfinsight$emorecogitem44[selfinsight$eraitem44==5] <- 1; selfinsight$emorecogitem44[selfinsight$eraitem44!=5] <- 0
selfinsight$emorecogitem45 <- NA; selfinsight$emorecogitem45[selfinsight$eraitem45==2] <- 1; selfinsight$emorecogitem45[selfinsight$eraitem45!=2] <- 0
selfinsight$emorecogitem46 <- NA; selfinsight$emorecogitem46[selfinsight$eraitem46==6] <- 1; selfinsight$emorecogitem46[selfinsight$eraitem46!=6] <- 0
selfinsight$emorecogitem47 <- NA; selfinsight$emorecogitem47[selfinsight$eraitem47==5] <- 1; selfinsight$emorecogitem47[selfinsight$eraitem47!=5] <- 0
selfinsight$emorecogitem48 <- NA; selfinsight$emorecogitem48[selfinsight$eraitem48==1] <- 1; selfinsight$emorecogitem48[selfinsight$eraitem48!=1] <- 0
selfinsight$emorecogitem49 <- NA; selfinsight$emorecogitem49[selfinsight$eraitem49==3] <- 1; selfinsight$emorecogitem49[selfinsight$eraitem49!=3] <- 0
selfinsight$emorecogitem50 <- NA; selfinsight$emorecogitem50[selfinsight$eraitem50==1] <- 1; selfinsight$emorecogitem50[selfinsight$eraitem50!=1] <- 0
selfinsight$emorecogitem51 <- NA; selfinsight$emorecogitem51[selfinsight$eraitem51==5] <- 1; selfinsight$emorecogitem51[selfinsight$eraitem51!=5] <- 0
selfinsight$emorecogitem52 <- NA; selfinsight$emorecogitem52[selfinsight$eraitem52==1] <- 1; selfinsight$emorecogitem52[selfinsight$eraitem52!=1] <- 0
selfinsight$emorecogitem53 <- NA; selfinsight$emorecogitem53[selfinsight$eraitem53==2] <- 1; selfinsight$emorecogitem53[selfinsight$eraitem53!=2] <- 0
selfinsight$emorecogitem54 <- NA; selfinsight$emorecogitem54[selfinsight$eraitem54==6] <- 1; selfinsight$emorecogitem54[selfinsight$eraitem54!=6] <- 0
selfinsight$emorecogitem55 <- NA; selfinsight$emorecogitem55[selfinsight$eraitem55==1] <- 1; selfinsight$emorecogitem55[selfinsight$eraitem55!=1] <- 0
selfinsight$emorecogitem56 <- NA; selfinsight$emorecogitem56[selfinsight$eraitem56==2] <- 1; selfinsight$emorecogitem56[selfinsight$eraitem56!=2] <- 0
selfinsight$emorecogitem57 <- NA; selfinsight$emorecogitem57[selfinsight$eraitem57==5] <- 1; selfinsight$emorecogitem57[selfinsight$eraitem57!=5] <- 0
selfinsight$emorecogitem58 <- NA; selfinsight$emorecogitem58[selfinsight$eraitem58==5] <- 1; selfinsight$emorecogitem58[selfinsight$eraitem58!=5] <- 0
selfinsight$emorecogitem59 <- NA; selfinsight$emorecogitem59[selfinsight$eraitem59==4] <- 1; selfinsight$emorecogitem59[selfinsight$eraitem59!=4] <- 0
selfinsight$emorecogitem60 <- NA; selfinsight$emorecogitem60[selfinsight$eraitem60==2] <- 1; selfinsight$emorecogitem60[selfinsight$eraitem60!=2] <- 0
selfinsight$emorecogitem61 <- NA; selfinsight$emorecogitem61[selfinsight$eraitem61==3] <- 1; selfinsight$emorecogitem61[selfinsight$eraitem61!=3] <- 0
selfinsight$emorecogitem62 <- NA; selfinsight$emorecogitem62[selfinsight$eraitem62==6] <- 1; selfinsight$emorecogitem62[selfinsight$eraitem62!=6] <- 0
selfinsight$emorecogitem63 <- NA; selfinsight$emorecogitem63[selfinsight$eraitem63==2] <- 1; selfinsight$emorecogitem63[selfinsight$eraitem63!=2] <- 0
selfinsight$emorecogitem64 <- NA; selfinsight$emorecogitem64[selfinsight$eraitem64==1] <- 1; selfinsight$emorecogitem64[selfinsight$eraitem64!=1] <- 0
selfinsight$emorecogitem65 <- NA; selfinsight$emorecogitem65[selfinsight$eraitem65==3] <- 1; selfinsight$emorecogitem65[selfinsight$eraitem65!=3] <- 0
selfinsight$emorecogitem66 <- NA; selfinsight$emorecogitem66[selfinsight$eraitem66==4] <- 1; selfinsight$emorecogitem66[selfinsight$eraitem66!=4] <- 0
selfinsight$emorecogitem67 <- NA; selfinsight$emorecogitem67[selfinsight$eraitem67==6] <- 1; selfinsight$emorecogitem67[selfinsight$eraitem67!=6] <- 0
selfinsight$emorecogitem68 <- NA; selfinsight$emorecogitem68[selfinsight$eraitem68==6] <- 1; selfinsight$emorecogitem68[selfinsight$eraitem68!=6] <- 0
selfinsight$emorecogitem69 <- NA; selfinsight$emorecogitem69[selfinsight$eraitem69==4] <- 1; selfinsight$emorecogitem69[selfinsight$eraitem69!=4] <- 0
selfinsight$emorecogitem70 <- NA; selfinsight$emorecogitem70[selfinsight$eraitem70==3] <- 1; selfinsight$emorecogitem70[selfinsight$eraitem70!=3] <- 0
selfinsight$emorecogitem71 <- NA; selfinsight$emorecogitem71[selfinsight$eraitem71==2] <- 1; selfinsight$emorecogitem71[selfinsight$eraitem71!=2] <- 0
selfinsight$emorecogitem72 <- NA; selfinsight$emorecogitem72[selfinsight$eraitem72==3] <- 1; selfinsight$emorecogitem72[selfinsight$eraitem72!=3] <- 0
selfinsight$ERAtotal <- with(selfinsight, (rowSums(cbind(emorecogitem1,emorecogitem2,emorecogitem3,emorecogitem4,emorecogitem5,emorecogitem6,emorecogitem7,emorecogitem8,emorecogitem9,emorecogitem10,emorecogitem11,emorecogitem12,emorecogitem13,emorecogitem14,emorecogitem15,emorecogitem16,emorecogitem17,emorecogitem18,emorecogitem19,emorecogitem20,emorecogitem21,emorecogitem22,emorecogitem23,emorecogitem24,emorecogitem25,emorecogitem26,emorecogitem27,emorecogitem28,emorecogitem29,emorecogitem30,emorecogitem31,emorecogitem32,emorecogitem33,emorecogitem34,emorecogitem35,emorecogitem36,emorecogitem37,emorecogitem38,emorecogitem39,emorecogitem40,emorecogitem41,emorecogitem42,emorecogitem43,emorecogitem44,emorecogitem45,emorecogitem46,emorecogitem47,emorecogitem48,emorecogitem49,emorecogitem50,emorecogitem51,emorecogitem52,emorecogitem53,emorecogitem54,emorecogitem55,emorecogitem56,emorecogitem57,emorecogitem58,emorecogitem59,emorecogitem60,emorecogitem61,emorecogitem62,emorecogitem63,emorecogitem64,emorecogitem65,emorecogitem66,emorecogitem67,emorecogitem68,emorecogitem69,emorecogitem70,emorecogitem71,emorecogitem72), na.rm=TRUE)))
table(selfinsight$ERAtotal)
selfinsight$ERApercent <- selfinsight$ERAtotal/72
#### calculate split half reliability of emotional abilities test
er.df <- with(selfinsight, data.frame(cbind(emorecogitem1,emorecogitem2,emorecogitem3,emorecogitem4,emorecogitem5,emorecogitem6,emorecogitem7,emorecogitem8,emorecogitem9,emorecogitem10,emorecogitem11,emorecogitem12,emorecogitem13,emorecogitem14,emorecogitem15,emorecogitem16,emorecogitem17,emorecogitem18,emorecogitem19,emorecogitem20,emorecogitem21,emorecogitem22,emorecogitem23,emorecogitem24,emorecogitem25,emorecogitem26,emorecogitem27,emorecogitem28,emorecogitem29,emorecogitem30,emorecogitem31,emorecogitem32,emorecogitem33,emorecogitem34,emorecogitem35,emorecogitem36,emorecogitem37,emorecogitem38,emorecogitem39,emorecogitem40,emorecogitem41,emorecogitem42,emorecogitem43,emorecogitem44,emorecogitem45,emorecogitem46,emorecogitem47,emorecogitem48,emorecogitem49,emorecogitem50,emorecogitem51,emorecogitem52,emorecogitem53,emorecogitem54,emorecogitem55,emorecogitem56,emorecogitem57,emorecogitem58,emorecogitem59,emorecogitem60,emorecogitem61,emorecogitem62,emorecogitem63,emorecogitem64,emorecogitem65,emorecogitem66,emorecogitem67,emorecogitem68,emorecogitem69,emorecogitem70,emorecogitem71,emorecogitem72)))
splithalf.r(er.df, sims = 1000, graph = TRUE)
#### if reliability is lower than .60, remove items with lowest item-total correlations until .60 is achieved (p.30 and p. 31 of pre-registration) 
#### note: reliability is .89 so there is no need to use this procedure

#### score the test of cognitive abilities
selfinsight$ravenitem1 <- NA; selfinsight$ravenitem1[selfinsight$Raven1==7] <- 1; selfinsight$ravenitem1[selfinsight$Raven1!=7] <- 0
selfinsight$ravenitem2 <- NA; selfinsight$ravenitem2[selfinsight$Raven2==4] <- 1; selfinsight$ravenitem2[selfinsight$Raven2!=4] <- 0
selfinsight$ravenitem3 <- NA; selfinsight$ravenitem3[selfinsight$Raven3==6] <- 1; selfinsight$ravenitem3[selfinsight$Raven3!=6] <- 0
selfinsight$ravenitem4 <- NA; selfinsight$ravenitem4[selfinsight$Raven4==2] <- 1; selfinsight$ravenitem4[selfinsight$Raven4!=2] <- 0
selfinsight$ravenitem5 <- NA; selfinsight$ravenitem5[selfinsight$Raven5==4] <- 1; selfinsight$ravenitem5[selfinsight$Raven5!=4] <- 0
selfinsight$ravenitem6 <- NA; selfinsight$ravenitem6[selfinsight$Raven6==7] <- 1; selfinsight$ravenitem6[selfinsight$Raven6!=7] <- 0
selfinsight$ravenitem7 <- NA; selfinsight$ravenitem7[selfinsight$Raven7==3] <- 1; selfinsight$ravenitem7[selfinsight$Raven7!=3] <- 0
selfinsight$ravenitem8 <- NA; selfinsight$ravenitem8[selfinsight$Raven8==8] <- 1; selfinsight$ravenitem8[selfinsight$Raven8!=8] <- 0
selfinsight$ravenitem9 <- NA; selfinsight$ravenitem9[selfinsight$Raven9==7] <- 1; selfinsight$ravenitem9[selfinsight$Raven9!=7] <- 0
selfinsight$ravenitem10 <- NA; selfinsight$ravenitem10[selfinsight$Raven10==3] <- 1; selfinsight$ravenitem10[selfinsight$Raven10!=3] <- 0
selfinsight$ravenitem11 <- NA; selfinsight$ravenitem11[selfinsight$Raven11==7] <- 1; selfinsight$ravenitem11[selfinsight$Raven11!=7] <- 0
selfinsight$ravenitem12 <- NA; selfinsight$ravenitem12[selfinsight$Raven12==5] <- 1; selfinsight$ravenitem12[selfinsight$Raven12!=5] <- 0
selfinsight$ravenitem13 <- NA; selfinsight$ravenitem13[selfinsight$Raven13==5] <- 1; selfinsight$ravenitem13[selfinsight$Raven13!=5] <- 0
selfinsight$ravenitem14 <- NA; selfinsight$ravenitem14[selfinsight$Raven14==4] <- 1; selfinsight$ravenitem14[selfinsight$Raven14!=4] <- 0
selfinsight$ravenitem15 <- NA; selfinsight$ravenitem15[selfinsight$Raven15==1] <- 1; selfinsight$ravenitem15[selfinsight$Raven15!=1] <- 0
selfinsight$raventotal <- with(selfinsight, (rowSums(cbind(ravenitem1,ravenitem2,ravenitem3,ravenitem4,ravenitem5,ravenitem6,ravenitem7,ravenitem8,ravenitem9,ravenitem10,ravenitem11,ravenitem12,ravenitem13,ravenitem14,ravenitem15), na.rm=TRUE)))
table(selfinsight$raventotal)
selfinsight$ravenpercent <- selfinsight$raventotal/15
#### calculate split half reliability of cognitive abilities test
#### first create a data frame with the 1,038 participants (out of 1,044) with valid cognitive ability data (as per the pre-registration) 
table(selfinsight$raventotalattempted)
selfinsight.valid.cognitive.ability <- subset(selfinsight, raventotalattempted>7 & ravenselfrating>-1)
raven.rel <- with(selfinsight.valid.cognitive.ability, data.frame(cbind(ravenitem1,ravenitem2,ravenitem3,ravenitem4,ravenitem5,ravenitem6,ravenitem7,ravenitem8,ravenitem9,ravenitem10,ravenitem11,ravenitem12,ravenitem13,ravenitem14,ravenitem15)))
splithalf.r(raven.rel, sims = 1000, graph = TRUE)
#### if reliability is lower than .60, remove items with lowest item-total correlations until .60 is achieved (p.30 and p. 31 of pre-registration) 
#### note: reliability is .79 so there is no need to use this procedure

#### EXAMINE EXISTENCE OF MATCHES AND MISMATCHES (pp. 30-31 of pre-registration) ####

### EMOTIONAL ABILITIES

### calculate standardized scores (z-scores) for each participant 
selfinsight$z.ERA <- (selfinsight$ERAtotal-36)/(sd(selfinsight$ERAtotal,na.rm=TRUE))
selfinsight$z.SV <- (selfinsight$ertselfrating-36)/(sd(selfinsight$ertselfrating,na.rm=TRUE))
with(selfinsight,describe(cbind(z.ERA, z.SV)))

### calculate the people who had accurate self-insight within half a SD
selfinsight$overestimater <- NA
selfinsight$overestimater[(selfinsight$z.ERA)<(selfinsight$z.SV-.5)] <- 1
selfinsight$overestimater[(selfinsight$z.ERA)>=(selfinsight$z.SV-.5)] <- 0
selfinsight$underestimater <- NA
selfinsight$underestimater[(selfinsight$z.ERA)>(selfinsight$z.SV+.5)] <- 1
selfinsight$underestimater[(selfinsight$z.ERA)<=(selfinsight$z.SV+.5)] <- 0
selfinsight$accurate <- NA; 
selfinsight$accurate[(selfinsight$z.ERA)>=(selfinsight$z.SV-.5) & (selfinsight$z.ERA)<=(selfinsight$z.SV+.5)] <- 1
selfinsight$accurate[(selfinsight$z.ERA)<(selfinsight$z.SV-.5) | (selfinsight$z.ERA)>(selfinsight$z.SV+.5)] <- 0

### display the people who were over- or under- or correct estimaters #
xtabs(~ overestimater, data = selfinsight)
xtabs(~ underestimater, data = selfinsight)
xtabs(~ accurate, data = selfinsight)

### COGNITIVE ABILITIES

### calculate standardized scores (z-scores) for each participant 
selfinsight.valid.cognitive.ability$z.raven <- (selfinsight.valid.cognitive.ability$raventotal-7.5)/(sd(selfinsight.valid.cognitive.ability$raventotal,na.rm=TRUE))
selfinsight.valid.cognitive.ability$z.ravenSV <- (selfinsight.valid.cognitive.ability$ravenselfrating-7.5)/(sd(selfinsight.valid.cognitive.ability$ravenselfrating,na.rm=TRUE))
with(selfinsight.valid.cognitive.ability,describe(cbind(z.raven, z.ravenSV)))

### calculate the people who had accurate self-insight within half a SD
selfinsight.valid.cognitive.ability$overestimater <- NA
selfinsight.valid.cognitive.ability$overestimater[(selfinsight.valid.cognitive.ability$z.raven)<(selfinsight.valid.cognitive.ability$z.ravenSV-.5)] <- 1
selfinsight.valid.cognitive.ability$overestimater[(selfinsight.valid.cognitive.ability$z.raven)>=(selfinsight.valid.cognitive.ability$z.ravenSV-.5)] <- 0
selfinsight.valid.cognitive.ability$underestimater <- NA
selfinsight.valid.cognitive.ability$underestimater[(selfinsight.valid.cognitive.ability$z.raven)>(selfinsight.valid.cognitive.ability$z.ravenSV+.5)] <- 1
selfinsight.valid.cognitive.ability$underestimater[(selfinsight.valid.cognitive.ability$z.raven)<=(selfinsight.valid.cognitive.ability$z.ravenSV+.5)] <- 0
selfinsight.valid.cognitive.ability$accurate <- NA; 
selfinsight.valid.cognitive.ability$accurate[(selfinsight.valid.cognitive.ability$z.raven)>=(selfinsight.valid.cognitive.ability$z.ravenSV-.5) & (selfinsight.valid.cognitive.ability$z.raven)<=(selfinsight.valid.cognitive.ability$z.ravenSV+.5)] <- 1
selfinsight.valid.cognitive.ability$accurate[(selfinsight.valid.cognitive.ability$z.raven)<(selfinsight.valid.cognitive.ability$z.ravenSV-.5) | (selfinsight.valid.cognitive.ability$z.raven)>(selfinsight.valid.cognitive.ability$z.ravenSV+.5)] <- 0

### display the people who were over- or under- or correct estimaters
xtabs(~ overestimater, data = selfinsight.valid.cognitive.ability)
xtabs(~ underestimater, data = selfinsight.valid.cognitive.ability)
xtabs(~ accurate, data = selfinsight.valid.cognitive.ability)

#### DESCRIPTIVES STATS AND (DESCRIPTIVE) CORRELATIONS ####

#### emotional abilities test scores and self-views of emotional abilities
with(selfinsight, describe(cbind(ERAtotal,ERAtotalattempted,ertselfrating)))

#### correlate emotional abilities test scores with self-views of emotional abilities
with(selfinsight, cor.test(ERAtotal,ertselfrating), use="pairwise", adjust="none")
with(selfinsight, corr.test(cbind(ERAtotal,ERAtotalattempted,ertselfrating)), use="pairwise", adjust="none")

#### cognitive abilities test scores and self-views of cognitive abilities
with(selfinsight.valid.cognitive.ability, describe(cbind(raventotal,raventotalattempted,ravenselfrating)))

#### correlate emotional and cognitive abilities test scores with self-views of emotional and cognitive abilities (note that this PAIRWISE and only includes observations with valid cognitive ability data as per the pre-registration)
with(selfinsight.valid.cognitive.ability, cor.test(raventotal,ravenselfrating), use="pairwise", adjust="none")
with(selfinsight.valid.cognitive.ability, corr.test(cbind(ERAtotal,ERAtotalattempted,ertselfrating,raventotal,raventotalattempted,ravenselfrating)), use="pairwise", adjust="none")

#### CALCULATE DESCRIPTIVES AND RELIABILITY OF ADJUSTMENT USING DATA MERGED DIFFERENTLY SO THAT EACH DIARY COMPLETED BY ANYONE IS A ROW ####

merged.data.long <- read.csv("Stage 2 - dairies long form.csv", header=TRUE, stringsAsFactors = FALSE)

#### calculate average duration across all diaries that were done
merged.data.long$Duration..in.seconds.<- as.numeric(merged.data.long$Duration..in.seconds.)
describe(merged.data.long$Duration..in.seconds.)
#### note: average duration of diary is 197.85 seconds (3.30 minutes)
#### note: but the max duration was 13536 seconds (225.60 minutes) indicating that the participant might have taken a break and left the survey open
meanduration <- mean(merged.data.long$Duration..in.seconds.); sdduration <- sd(merged.data.long$Duration..in.seconds.); cutoffduration <- meanduration+(3.29*sdduration)
meanduration; sdduration; cutoffduration
merged.data.long.nodurationoutlier <- subset(merged.data.long, Duration..in.seconds.<cutoffduration)
describe(merged.data.long.nodurationoutlier$Duration..in.seconds.)
#### average duration without outliers is 147.34 seconds (2.46 minutes)

#### calculate descriptives and reliability of adjustment measures across all diaries that were done by all participants

merged.data.long$swls_swls_1 <- as.numeric(merged.data.long$swls_swls_1); merged.data.long$swls_swls_2 <- as.numeric(merged.data.long$swls_swls_2); merged.data.long$swls_swls_3 <- as.numeric(merged.data.long$swls_swls_3); merged.data.long$swls_swls_4 <- as.numeric(merged.data.long$swls_swls_4); merged.data.long$swls_swls_5 <- as.numeric(merged.data.long$swls_swls_5)
merged.data.long$life.satisfaction <- with(merged.data.long,(rowMeans(cbind(swls_swls_1,swls_swls_2,swls_swls_3,swls_swls_4,swls_swls_5),na.rm=TRUE)))
with(merged.data.long, (psych::alpha(data.frame(swls_swls_1,swls_swls_2,swls_swls_3,swls_swls_4,swls_swls_5),na.rm=TRUE)))

merged.data.long$institutional_institutional_1 <- as.numeric(merged.data.long$institutional_institutional_1); merged.data.long$institutional_institutional_2 <- as.numeric(merged.data.long$institutional_institutional_2); merged.data.long$institutional_institutional_3 <- as.numeric(merged.data.long$institutional_institutional_3); merged.data.long$institutional_institutional_4 <- as.numeric(merged.data.long$institutional_institutional_4); merged.data.long$institutional_institutional_5 <- as.numeric(merged.data.long$institutional_institutional_5)
merged.data.long$career.satisfaction <- with(merged.data.long,(rowMeans(cbind(institutional_institutional_1,institutional_institutional_2,institutional_institutional_3,institutional_institutional_4,institutional_institutional_5),na.rm=TRUE)))
with(merged.data.long, (psych::alpha(data.frame(institutional_institutional_1,institutional_institutional_2,institutional_institutional_3,institutional_institutional_4,institutional_institutional_5),na.rm=TRUE)))

merged.data.long$interpersonal_interpersonal_1 <- as.numeric(merged.data.long$interpersonal_interpersonal_1); merged.data.long$interpersonal_interpersonal_2 <- as.numeric(merged.data.long$interpersonal_interpersonal_2); merged.data.long$interpersonal_interpersonal_3 <- as.numeric(merged.data.long$interpersonal_interpersonal_3); merged.data.long$interpersonal_interpersonal_4 <- as.numeric(merged.data.long$interpersonal_interpersonal_4)
merged.data.long$r.interpersonal_interpersonal_3 <- 6-merged.data.long$interpersonal_interpersonal_3; merged.data.long$r.interpersonal_interpersonal_4 <- 6-merged.data.long$interpersonal_interpersonal_4 
merged.data.long$relationship.satisfaction <- with(merged.data.long,(rowMeans(cbind(interpersonal_interpersonal_1,interpersonal_interpersonal_2,r.interpersonal_interpersonal_3,r.interpersonal_interpersonal_4),na.rm=TRUE)))
with(merged.data.long, (psych::alpha(data.frame(interpersonal_interpersonal_1,interpersonal_interpersonal_2,r.interpersonal_interpersonal_3,r.interpersonal_interpersonal_4),na.rm=TRUE)))

with(merged.data.long, describe(cbind(life.satisfaction,career.satisfaction,relationship.satisfaction)))

#### examine missing data in the diaries
missing.lifesatisfactionitems <- subset(merged.data.long, select=swls_swls_1:swls_swls_5); sum(is.na(missing.lifesatisfactionitems))
missing.careersatisfactionitems <- subset(merged.data.long, select=institutional_institutional_1:institutional_institutional_5); sum(is.na(missing.careersatisfactionitems))
missing.relationshipsatisfactionitems <- subset(merged.data.long, select=interpersonal_interpersonal_1:interpersonal_interpersonal_4); sum(is.na(missing.relationshipsatisfactionitems))
missing.lifesatisfactioncomposite <- subset(merged.data.long, select=life.satisfaction); sum(is.na(missing.lifesatisfactioncomposite))
missing.careersatisfactioncomposite <- subset(merged.data.long, select=career.satisfaction); sum(is.na(missing.careersatisfactioncomposite))
missing.relationshipsatisfactioncomposite <- subset(merged.data.long, select=relationship.satisfaction); sum(is.na(missing.relationshipsatisfactioncomposite))
merged.data.long.missing.career.composite <- subset(merged.data.long, is.na(career.satisfaction)) 
table(merged.data.long.missing.career.composite$life.satisfaction)
table(merged.data.long.missing.career.composite$career.satisfaction)
table(merged.data.long.missing.career.composite$relationship.satisfaction)

#### ANALYSES WITH ERA (EMOTION RECOGNITION ABILITY) ####

#### CENTER VARIABLES

selfinsight$c.ERAtotal <- selfinsight$ERAtotal-36; selfinsight$r.c.ERAtotal <- 0-selfinsight$c.ERAtotal
selfinsight$c.ertselfrating <- selfinsight$ertselfrating-36; selfinsight$r.c.ertselfrating <- 0-selfinsight$c.ertselfrating
with(selfinsight,describe(cbind(c.ERAtotal,r.c.ERAtotal,ERAtotal,c.ertselfrating,r.c.ertselfrating,ertselfrating)))

#### ERA AND LIFE SATISFACTION ####

# calculate life satisfaction score for each diary

selfinsight$swls_swls_1_1 <- as.numeric(selfinsight$swls_swls_1_1); selfinsight$swls_swls_2_1 <- as.numeric(selfinsight$swls_swls_2_1); selfinsight$swls_swls_3_1 <- as.numeric(selfinsight$swls_swls_3_1); selfinsight$swls_swls_4_1 <- as.numeric(selfinsight$swls_swls_4_1); selfinsight$swls_swls_5_1 <- as.numeric(selfinsight$swls_swls_5_1)
selfinsight$swls.diary1 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_1,swls_swls_2_1,swls_swls_3_1,swls_swls_4_1,swls_swls_5_1),na.rm=TRUE)))
selfinsight$swls_swls_1_2 <- as.numeric(selfinsight$swls_swls_1_2); selfinsight$swls_swls_2_2 <- as.numeric(selfinsight$swls_swls_2_2); selfinsight$swls_swls_3_2 <- as.numeric(selfinsight$swls_swls_3_2); selfinsight$swls_swls_4_2 <- as.numeric(selfinsight$swls_swls_4_2); selfinsight$swls_swls_5_2 <- as.numeric(selfinsight$swls_swls_5_2)
selfinsight$swls.diary2 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_2,swls_swls_2_2,swls_swls_3_2,swls_swls_4_2,swls_swls_5_2),na.rm=TRUE)))
selfinsight$swls_swls_1_3 <- as.numeric(selfinsight$swls_swls_1_3); selfinsight$swls_swls_2_3 <- as.numeric(selfinsight$swls_swls_2_3); selfinsight$swls_swls_3_3 <- as.numeric(selfinsight$swls_swls_3_3); selfinsight$swls_swls_4_3 <- as.numeric(selfinsight$swls_swls_4_3); selfinsight$swls_swls_5_3 <- as.numeric(selfinsight$swls_swls_5_3)
selfinsight$swls.diary3 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_3,swls_swls_2_3,swls_swls_3_3,swls_swls_4_3,swls_swls_5_3),na.rm=TRUE)))
selfinsight$swls_swls_1_4 <- as.numeric(selfinsight$swls_swls_1_4); selfinsight$swls_swls_2_4 <- as.numeric(selfinsight$swls_swls_2_4); selfinsight$swls_swls_3_4 <- as.numeric(selfinsight$swls_swls_3_4); selfinsight$swls_swls_4_4 <- as.numeric(selfinsight$swls_swls_4_4); selfinsight$swls_swls_5_4 <- as.numeric(selfinsight$swls_swls_5_4)
selfinsight$swls.diary4 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_4,swls_swls_2_4,swls_swls_3_4,swls_swls_4_4,swls_swls_5_4),na.rm=TRUE)))
selfinsight$swls_swls_1_5 <- as.numeric(selfinsight$swls_swls_1_5); selfinsight$swls_swls_2_5 <- as.numeric(selfinsight$swls_swls_2_5); selfinsight$swls_swls_3_5 <- as.numeric(selfinsight$swls_swls_3_5); selfinsight$swls_swls_4_5 <- as.numeric(selfinsight$swls_swls_4_5); selfinsight$swls_swls_5_5 <- as.numeric(selfinsight$swls_swls_5_5)
selfinsight$swls.diary5 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_5,swls_swls_2_5,swls_swls_3_5,swls_swls_4_5,swls_swls_5_5),na.rm=TRUE)))
selfinsight$swls_swls_1_6 <- as.numeric(selfinsight$swls_swls_1_6); selfinsight$swls_swls_2_6 <- as.numeric(selfinsight$swls_swls_2_6); selfinsight$swls_swls_3_6 <- as.numeric(selfinsight$swls_swls_3_6); selfinsight$swls_swls_4_6 <- as.numeric(selfinsight$swls_swls_4_6); selfinsight$swls_swls_5_6 <- as.numeric(selfinsight$swls_swls_5_6)
selfinsight$swls.diary6 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_6,swls_swls_2_6,swls_swls_3_6,swls_swls_4_6,swls_swls_5_6),na.rm=TRUE)))
selfinsight$swls_swls_1_7 <- as.numeric(selfinsight$swls_swls_1_7); selfinsight$swls_swls_2_7 <- as.numeric(selfinsight$swls_swls_2_7); selfinsight$swls_swls_3_7 <- as.numeric(selfinsight$swls_swls_3_7); selfinsight$swls_swls_4_7 <- as.numeric(selfinsight$swls_swls_4_7); selfinsight$swls_swls_5_7 <- as.numeric(selfinsight$swls_swls_5_7)
selfinsight$swls.diary7 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_7,swls_swls_2_7,swls_swls_3_7,swls_swls_4_7,swls_swls_5_7),na.rm=TRUE)))
selfinsight$swls_swls_1_8 <- as.numeric(selfinsight$swls_swls_1_8); selfinsight$swls_swls_2_8 <- as.numeric(selfinsight$swls_swls_2_8); selfinsight$swls_swls_3_8 <- as.numeric(selfinsight$swls_swls_3_8); selfinsight$swls_swls_4_8 <- as.numeric(selfinsight$swls_swls_4_8); selfinsight$swls_swls_5_8 <- as.numeric(selfinsight$swls_swls_5_8)
selfinsight$swls.diary8 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_8,swls_swls_2_8,swls_swls_3_8,swls_swls_4_8,swls_swls_5_8),na.rm=TRUE)))
selfinsight$swls_swls_1_9 <- as.numeric(selfinsight$swls_swls_1_9); selfinsight$swls_swls_2_9 <- as.numeric(selfinsight$swls_swls_2_9); selfinsight$swls_swls_3_9 <- as.numeric(selfinsight$swls_swls_3_9); selfinsight$swls_swls_4_9 <- as.numeric(selfinsight$swls_swls_4_9); selfinsight$swls_swls_5_9 <- as.numeric(selfinsight$swls_swls_5_9)
selfinsight$swls.diary9 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_9,swls_swls_2_9,swls_swls_3_9,swls_swls_4_9,swls_swls_5_9),na.rm=TRUE)))
selfinsight$swls_swls_1_0 <- as.numeric(selfinsight$swls_swls_1_0); selfinsight$swls_swls_2_0 <- as.numeric(selfinsight$swls_swls_2_0); selfinsight$swls_swls_3_0 <- as.numeric(selfinsight$swls_swls_3_0); selfinsight$swls_swls_4_0 <- as.numeric(selfinsight$swls_swls_4_0); selfinsight$swls_swls_5_0 <- as.numeric(selfinsight$swls_swls_5_0)
selfinsight$swls.diary0 <- with(selfinsight,(rowMeans(cbind(swls_swls_1_0,swls_swls_2_0,swls_swls_3_0,swls_swls_4_0,swls_swls_5_0),na.rm=TRUE)))

# average scores for all diaries
selfinsight$life.satisfaction <- with(selfinsight,(rowMeans(cbind(swls.diary1,swls.diary2,swls.diary3,swls.diary4,swls.diary5,swls.diary6,swls.diary7,swls.diary8,swls.diary9,swls.diary0),na.rm=TRUE)))
with(selfinsight, describe(cbind(life.satisfaction,swls.diary1,swls.diary2,swls.diary3,swls.diary4,swls.diary5,swls.diary6,swls.diary7,swls.diary8,swls.diary9,swls.diary0)))

#### REMOVE OUTLIERS (p. 32 of pre-registration)

#### UNIVARIATE OUTLIERS
#### observations 3.29 standard deviations above and below the mean (Tabachnick & Fidell)

meanERAtotal <- mean(selfinsight$ERAtotal); sdERAtotal <- SD(selfinsight$ERAtotal); lowcutoffERAtotal <- meanERAtotal-(3.29*sdERAtotal); highcutoffERAtotal <- meanERAtotal+(3.29*sdERAtotal)
table(selfinsight$ERAtotal)
meanERAtotal; sdERAtotal; lowcutoffERAtotal; highcutoffERAtotal 
#### note: 11 observations are LOW outliers on the emotion recognition test
#### as per the pre-registration, these data will be removed from the analyses
meanertselfrating <- mean(selfinsight$ertselfrating); sdertselfrating <- SD(selfinsight$ertselfrating); lowcutoffertselfrating <- meanertselfrating-(3.29*sdertselfrating); highcutoffertselfrating <- meanertselfrating+(3.29*sdertselfrating)
table(selfinsight$ertselfrating)
meanertselfrating; sdertselfrating; lowcutoffertselfrating; highcutoffertselfrating 
#### note: 15 observations are LOW outliers on self-views about emotion recognition
#### as per the pre-registration, these analyses will be removed from the analyses
meanlife.satisfaction <- mean(selfinsight$life.satisfaction); sdlife.satisfaction <- SD(selfinsight$life.satisfaction); lowcutofflife.satisfaction <- meanlife.satisfaction-(3.29*sdlife.satisfaction); highcutofflife.satisfaction <- meanlife.satisfaction+(3.29*sdlife.satisfaction)
meanlife.satisfaction; sdlife.satisfaction; lowcutofflife.satisfaction; highcutofflife.satisfaction 
#### note: NO observations are outliers on life satisfaction (cut offs are actually outside the possible range for the measure)
selfinsight.emo.lifesat.temp1 <- subset(selfinsight, ERAtotal<highcutoffERAtotal)
selfinsight.emo.lifesat.temp2 <- subset(selfinsight.emo.lifesat.temp1, ERAtotal>lowcutoffERAtotal)
selfinsight.emo.lifesat.temp3 <- subset(selfinsight.emo.lifesat.temp2, ertselfrating<highcutoffertselfrating)
selfinsight.emo.lifesat.temp4 <- subset(selfinsight.emo.lifesat.temp3, ertselfrating>lowcutoffertselfrating)
selfinsight.emo.lifesat.temp5 <- subset(selfinsight.emo.lifesat.temp4, life.satisfaction<highcutofflife.satisfaction)
selfinsight.emo.lifesat <- subset(selfinsight.emo.lifesat.temp5, life.satisfaction>lowcutofflife.satisfaction)

#### multivariate outliers
#### observations with values that exceed the cut-offs recommended by Bollen and Jackman98 and Belsley, Kuh, and Welsch on three indices: Cook's D, Hat values, and difference in fits (DFFITS)
#### leverage measures
selfinsight.emo.lifesat$out <- FALSE
selfinsight.emo.lifesat$c.ERAtotal2 <- (selfinsight.emo.lifesat$c.ERAtotal)^2
selfinsight.emo.lifesat$c.ertselfrating2 <- (selfinsight.emo.lifesat$c.ertselfrating)^2
selfinsight.emo.lifesat$c.ERAtotal_c.ertselfrating <- selfinsight.emo.lifesat$c.ERAtotal*selfinsight.emo.lifesat$c.ertselfrating
full_model <- lm(life.satisfaction~c.ERAtotal+c.ertselfrating+c.ERAtotal2+c.ERAtotal_c.ertselfrating+c.ertselfrating2, data=selfinsight.emo.lifesat, na.action=na.exclude); summary(full_model)
infl <- influence.measures(full_model)
# function to get leverage measures (Bollen & Jackman, 1985)
selfinsight.emo.lifesat$out <- apply(infl$is.inf[, c("dffit", "cook.d", "hat")], 1, sum) == 3
# if a case has significant cooks d, dffit, AND hat, then it gets flagged as an outlier
n.out <- sum(na.omit(selfinsight.emo.lifesat$out) == TRUE)
# count the number of outliers identified by the previous line
n.out
# no outliers in data
# if outliers existed, they would be taken out (via data=accuracy[selfinsight.emo.lifesat.temp6$out==FALSE, ])

#### multicollinearity
ols_vif_tol(full_model)

#### testing assumptions of OLS regression ####

# normality of residuals
ols_plot_resid_qq(full_model)
ols_test_normality(full_model)
ols_plot_resid_hist(full_model)

# heteroskedasticity
bptest(full_model)

#### the data violate assumptions of homoscedasticity; use robust standard errors

#### correlations between variables
with(selfinsight.emo.lifesat, corr.test(cbind(c.ERAtotal,c.ertselfrating,life.satisfaction)))

#### test of significance of polynomial terms
summary(full_model.1.a <- lm(life.satisfaction~c.ERAtotal+c.ertselfrating, data=selfinsight.emo.lifesat, na.action=na.exclude))
summary(full_model.1.b <- lm(life.satisfaction~c.ERAtotal+c.ertselfrating+c.ERAtotal2+c.ERAtotal_c.ertselfrating+c.ertselfrating2, data=selfinsight.emo.lifesat, na.action=na.exclude))
anova(full_model.1.a, full_model.1.b)

#### test for self-insight with polynomial regression analyses and RSA
RSA.emotion.lifesat<-RSA(life.satisfaction~c.ERAtotal*c.ertselfrating,data=selfinsight.emo.lifesat,
                center=FALSE,
                scale=FALSE,
                na.rm=TRUE,
                out.rm=NULL,
                breakline=FALSE, 
                models="default",
                cubic=FALSE,
                verbose=TRUE,
                estimator="MLR",
                se="robust",
                missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.emotion.lifesat)
getPar(RSA.emotion.lifesat)
# plot RSA graph 
plot(RSA.emotion.lifesat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1", "PA2"), 
     project = c("PA1", "PA2", "LOC"),        
     xlim    = c(-36,36),   
     ylim    = c(-36,36),    
     xlab = "Emotional Ability", 
     ylab = "Self View", 
     zlab = "Life Satisfaction"
)
# plot differently to see it better
RSA.emotion.lifesat.reverse<-RSA(life.satisfaction~r.c.ERAtotal*r.c.ertselfrating,data=selfinsight.emo.lifesat,
                         center=FALSE,
                         scale=FALSE,
                         na.rm=TRUE,
                         out.rm=NULL,
                         breakline=FALSE, 
                         models="default",
                         cubic=FALSE,
                         verbose=TRUE,
                         estimator="MLR",
                         se="robust",
                         missing='listwise'
)
# plot reverse RSA graph 
plot(RSA.emotion.lifesat.reverse,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-36,36),   
     ylim    = c(-36,36),    
     xlab = "Emotional Ability (reversed)", 
     ylab = "Self View (reversed)", 
     zlab = "Life Satisfaction"
)

#### ROBUSTNESS CHECK WITH COMMON GRAND MEAN CENTERING
commonmean <- (mean(selfinsight.emo.lifesat$ERAtotal)+mean(selfinsight.emo.lifesat$ertselfrating))/2
selfinsight.emo.lifesat$cc.ERAtotal <- (selfinsight.emo.lifesat$ERAtotal)-commonmean
selfinsight.emo.lifesat$cc.ertselfrating <- (selfinsight.emo.lifesat$ertselfrating)-commonmean
with(selfinsight.emo.lifesat,describe(cbind(cc.ERAtotal,ERAtotal,cc.ertselfrating,ertselfrating,commonmean)))
#### test for self-insight with polynomial regression analyses and RSA
RSA.emotion.lifesat.cc<-RSA(life.satisfaction~cc.ERAtotal*cc.ertselfrating,data=selfinsight.emo.lifesat,
                  center=FALSE,
                  scale=FALSE,
                  na.rm=TRUE,
                  out.rm=NULL,
                  breakline=FALSE, 
                  models="default",
                  cubic=FALSE,
                  verbose=TRUE,
                  estimator="MLR",
                  se="robust",
                  missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.emotion.lifesat.cc)
getPar(RSA.emotion.lifesat.cc)
# plot RSA graph 
plot(RSA.emotion.lifesat.cc,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-52.47,19.53),   
     ylim    = c(-52.47,19.53),    
     xlab = "Emotional Ability", 
     ylab = "Self View", 
     zlab = "Life Satisfaction"
)

#### ERA AND CAREER SATISFACTION ####

# calculate career satisfaction score for each diary

selfinsight$institutional_institutional_1_1 <- as.numeric(selfinsight$institutional_institutional_1_1); selfinsight$institutional_institutional_2_1 <- as.numeric(selfinsight$institutional_institutional_2_1); selfinsight$institutional_institutional_3_1 <- as.numeric(selfinsight$institutional_institutional_3_1); selfinsight$institutional_institutional_4_1 <- as.numeric(selfinsight$institutional_institutional_4_1); selfinsight$institutional_institutional_5_1 <- as.numeric(selfinsight$institutional_institutional_5_1)
selfinsight$institutional.diary1 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_1,institutional_institutional_2_1,institutional_institutional_3_1,institutional_institutional_4_1,institutional_institutional_5_1),na.rm=TRUE)))
selfinsight$institutional_institutional_1_2 <- as.numeric(selfinsight$institutional_institutional_1_2); selfinsight$institutional_institutional_2_2 <- as.numeric(selfinsight$institutional_institutional_2_2); selfinsight$institutional_institutional_3_2 <- as.numeric(selfinsight$institutional_institutional_3_2); selfinsight$institutional_institutional_4_2 <- as.numeric(selfinsight$institutional_institutional_4_2); selfinsight$institutional_institutional_5_2 <- as.numeric(selfinsight$institutional_institutional_5_2)
selfinsight$institutional.diary2 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_2,institutional_institutional_2_2,institutional_institutional_3_2,institutional_institutional_4_2,institutional_institutional_5_2),na.rm=TRUE)))
selfinsight$institutional_institutional_1_3 <- as.numeric(selfinsight$institutional_institutional_1_3); selfinsight$institutional_institutional_2_3 <- as.numeric(selfinsight$institutional_institutional_2_3); selfinsight$institutional_institutional_3_3 <- as.numeric(selfinsight$institutional_institutional_3_3); selfinsight$institutional_institutional_4_3 <- as.numeric(selfinsight$institutional_institutional_4_3); selfinsight$institutional_institutional_5_3 <- as.numeric(selfinsight$institutional_institutional_5_3)
selfinsight$institutional.diary3 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_3,institutional_institutional_2_3,institutional_institutional_3_3,institutional_institutional_4_3,institutional_institutional_5_3),na.rm=TRUE)))
selfinsight$institutional_institutional_1_4 <- as.numeric(selfinsight$institutional_institutional_1_4); selfinsight$institutional_institutional_2_4 <- as.numeric(selfinsight$institutional_institutional_2_4); selfinsight$institutional_institutional_3_4 <- as.numeric(selfinsight$institutional_institutional_3_4); selfinsight$institutional_institutional_4_4 <- as.numeric(selfinsight$institutional_institutional_4_4); selfinsight$institutional_institutional_5_4 <- as.numeric(selfinsight$institutional_institutional_5_4)
selfinsight$institutional.diary4 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_4,institutional_institutional_2_4,institutional_institutional_3_4,institutional_institutional_4_4,institutional_institutional_5_4),na.rm=TRUE)))
selfinsight$institutional_institutional_1_5 <- as.numeric(selfinsight$institutional_institutional_1_5); selfinsight$institutional_institutional_2_5 <- as.numeric(selfinsight$institutional_institutional_2_5); selfinsight$institutional_institutional_3_5 <- as.numeric(selfinsight$institutional_institutional_3_5); selfinsight$institutional_institutional_4_5 <- as.numeric(selfinsight$institutional_institutional_4_5); selfinsight$institutional_institutional_5_5 <- as.numeric(selfinsight$institutional_institutional_5_5)
selfinsight$institutional.diary5 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_5,institutional_institutional_2_5,institutional_institutional_3_5,institutional_institutional_4_5,institutional_institutional_5_5),na.rm=TRUE)))
selfinsight$institutional_institutional_1_6 <- as.numeric(selfinsight$institutional_institutional_1_6); selfinsight$institutional_institutional_2_6 <- as.numeric(selfinsight$institutional_institutional_2_6); selfinsight$institutional_institutional_3_6 <- as.numeric(selfinsight$institutional_institutional_3_6); selfinsight$institutional_institutional_4_6 <- as.numeric(selfinsight$institutional_institutional_4_6); selfinsight$institutional_institutional_5_6 <- as.numeric(selfinsight$institutional_institutional_5_6)
selfinsight$institutional.diary6 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_6,institutional_institutional_2_6,institutional_institutional_3_6,institutional_institutional_4_6,institutional_institutional_5_6),na.rm=TRUE)))
selfinsight$institutional_institutional_1_7 <- as.numeric(selfinsight$institutional_institutional_1_7); selfinsight$institutional_institutional_2_7 <- as.numeric(selfinsight$institutional_institutional_2_7); selfinsight$institutional_institutional_3_7 <- as.numeric(selfinsight$institutional_institutional_3_7); selfinsight$institutional_institutional_4_7 <- as.numeric(selfinsight$institutional_institutional_4_7); selfinsight$institutional_institutional_5_7 <- as.numeric(selfinsight$institutional_institutional_5_7)
selfinsight$institutional.diary7 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_7,institutional_institutional_2_7,institutional_institutional_3_7,institutional_institutional_4_7,institutional_institutional_5_7),na.rm=TRUE)))
selfinsight$institutional_institutional_1_8 <- as.numeric(selfinsight$institutional_institutional_1_8); selfinsight$institutional_institutional_2_8 <- as.numeric(selfinsight$institutional_institutional_2_8); selfinsight$institutional_institutional_3_8 <- as.numeric(selfinsight$institutional_institutional_3_8); selfinsight$institutional_institutional_4_8 <- as.numeric(selfinsight$institutional_institutional_4_8); selfinsight$institutional_institutional_5_8 <- as.numeric(selfinsight$institutional_institutional_5_8)
selfinsight$institutional.diary8 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_8,institutional_institutional_2_8,institutional_institutional_3_8,institutional_institutional_4_8,institutional_institutional_5_8),na.rm=TRUE)))
selfinsight$institutional_institutional_1_9 <- as.numeric(selfinsight$institutional_institutional_1_9); selfinsight$institutional_institutional_2_9 <- as.numeric(selfinsight$institutional_institutional_2_9); selfinsight$institutional_institutional_3_9 <- as.numeric(selfinsight$institutional_institutional_3_9); selfinsight$institutional_institutional_4_9 <- as.numeric(selfinsight$institutional_institutional_4_9); selfinsight$institutional_institutional_5_9 <- as.numeric(selfinsight$institutional_institutional_5_9)
selfinsight$institutional.diary9 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_9,institutional_institutional_2_9,institutional_institutional_3_9,institutional_institutional_4_9,institutional_institutional_5_9),na.rm=TRUE)))
selfinsight$institutional_institutional_1_0 <- as.numeric(selfinsight$institutional_institutional_1_0); selfinsight$institutional_institutional_2_0 <- as.numeric(selfinsight$institutional_institutional_2_0); selfinsight$institutional_institutional_3_0 <- as.numeric(selfinsight$institutional_institutional_3_0); selfinsight$institutional_institutional_4_0 <- as.numeric(selfinsight$institutional_institutional_4_0); selfinsight$institutional_institutional_5_0 <- as.numeric(selfinsight$institutional_institutional_5_0)
selfinsight$institutional.diary0 <- with(selfinsight,(rowMeans(cbind(institutional_institutional_1_0,institutional_institutional_2_0,institutional_institutional_3_0,institutional_institutional_4_0,institutional_institutional_5_0),na.rm=TRUE)))

# average scores for all diaries
selfinsight$career.satisfaction <- with(selfinsight,(rowMeans(cbind(institutional.diary1,institutional.diary2,institutional.diary3,institutional.diary4,institutional.diary5,institutional.diary6,institutional.diary7,institutional.diary8,institutional.diary9,institutional.diary0),na.rm=TRUE)))
with(selfinsight, describe(cbind(career.satisfaction,institutional.diary1,institutional.diary2,institutional.diary3,institutional.diary4,institutional.diary5,institutional.diary6,institutional.diary7,institutional.diary8,institutional.diary9,institutional.diary0)))

#### REMOVE OUTLIERS (p. 32 of pre-registration)

#### UNIVARIATE OUTLIERS
#### observations 3.29 standard deviations above and below the mean (Tabachnick & Fidell)

meanERAtotal <- mean(selfinsight$ERAtotal); sdERAtotal <- SD(selfinsight$ERAtotal); lowcutoffERAtotal <- meanERAtotal-(3.29*sdERAtotal); highcutoffERAtotal <- meanERAtotal+(3.29*sdERAtotal)
table(selfinsight$ERAtotal)
meanERAtotal; sdERAtotal; lowcutoffERAtotal; highcutoffERAtotal 
#### note: 11 observations are LOW outliers on the emotion recognition test
#### as per the pre-registration, these analyses will be removed from the analyses
meanertselfrating <- mean(selfinsight$ertselfrating); sdertselfrating <- SD(selfinsight$ertselfrating); lowcutoffertselfrating <- meanertselfrating-(3.29*sdertselfrating); highcutoffertselfrating <- meanertselfrating+(3.29*sdertselfrating)
table(selfinsight$ertselfrating)
meanertselfrating; sdertselfrating; lowcutoffertselfrating; highcutoffertselfrating 
#### note: 15 observations are LOW outliers on self-views about emotion recognition
#### as per the pre-registration, these analyses will be removed from the analyses
meancareer.satisfaction <- mean(selfinsight$career.satisfaction); sdcareer.satisfaction <- SD(selfinsight$career.satisfaction); lowcutoffcareer.satisfaction <- meancareer.satisfaction-(3.29*sdcareer.satisfaction); highcutoffcareer.satisfaction <- meancareer.satisfaction+(3.29*sdcareer.satisfaction)
meancareer.satisfaction; sdcareer.satisfaction; lowcutoffcareer.satisfaction; highcutoffcareer.satisfaction 
#### note: NO observations are outliers on career satisfaction (cut offs are actually outside the possible range for the measure)
selfinsight.emo.careersat.temp1 <- subset(selfinsight, ERAtotal<highcutoffERAtotal)
selfinsight.emo.careersat.temp2 <- subset(selfinsight.emo.careersat.temp1, ERAtotal>lowcutoffERAtotal)
selfinsight.emo.careersat.temp3 <- subset(selfinsight.emo.careersat.temp2, ertselfrating<highcutoffertselfrating)
selfinsight.emo.careersat.temp4 <- subset(selfinsight.emo.careersat.temp3, ertselfrating>lowcutoffertselfrating)
selfinsight.emo.careersat.temp5 <- subset(selfinsight.emo.careersat.temp4, career.satisfaction<highcutoffcareer.satisfaction)
selfinsight.emo.careersat <- subset(selfinsight.emo.careersat.temp5, career.satisfaction>lowcutoffcareer.satisfaction)

#### multivariate outliers
#### observations with values that exceed the cut-offs recommended by Bollen and Jackman98 and Belsley, Kuh, and Welsch on three indices: Cook's D, Hat values, and difference in fits (DFFITS)
#### leverage measures
selfinsight.emo.careersat$out <- FALSE
selfinsight.emo.careersat$c.ERAtotal2 <- (selfinsight.emo.careersat$c.ERAtotal)^2
selfinsight.emo.careersat$c.ertselfrating2 <- (selfinsight.emo.careersat$c.ertselfrating)^2
selfinsight.emo.careersat$c.ERAtotal_c.ertselfrating <- selfinsight.emo.careersat$c.ERAtotal*selfinsight.emo.careersat$c.ertselfrating
full_model <- lm(career.satisfaction~c.ERAtotal+c.ertselfrating+c.ERAtotal2+c.ERAtotal_c.ertselfrating+c.ertselfrating2, data=selfinsight.emo.careersat, na.action=na.exclude); summary(full_model)
infl <- influence.measures(full_model)
# function to get leverage measures (Bollen & Jackman, 1985)
selfinsight.emo.careersat$out <- apply(infl$is.inf[, c("dffit", "cook.d", "hat")], 1, sum) == 3
# if a case has significant cooks d, dffit, AND hat, then it gets flagged as an outlier
n.out <- sum(na.omit(selfinsight.emo.careersat$out) == TRUE)
# count the number of outliers identified by the previous line
n.out
# no outliers in data
# if outliers existed, they would be taken out (via data=accuracy[selfinsight.emo.careersat.temp6$out==FALSE, ])

#### multicollinearity
ols_vif_tol(full_model)

#### testing assumptions of OLS regression ####

# normality of residuals
ols_plot_resid_qq(full_model)
ols_test_normality(full_model)
ols_plot_resid_hist(full_model)

# heteroskedasticity
bptest(full_model)

#### the data violate assumptions of homoscedasticity; use robust standard errors

#### test of significance of polynomial terms
summary(full_model.2.a <- lm(career.satisfaction~c.ERAtotal+c.ertselfrating, data=selfinsight.emo.careersat, na.action=na.exclude))
summary(full_model.2.b <- lm(career.satisfaction~c.ERAtotal+c.ertselfrating+c.ERAtotal2+c.ERAtotal_c.ertselfrating+c.ertselfrating2, data=selfinsight.emo.careersat, na.action=na.exclude))
anova(full_model.2.a, full_model.2.b)

#### test for self-insight with polynomial regression analyses and RSA
RSA.emotion.careersat<-RSA(career.satisfaction~c.ERAtotal*c.ertselfrating,data=selfinsight.emo.careersat,
                         center=FALSE,
                         scale=FALSE,
                         na.rm=TRUE,
                         out.rm=NULL,
                         breakline=FALSE, 
                         models="default",
                         cubic=FALSE,
                         verbose=TRUE,
                         estimator="MLR",
                         se="robust",
                         missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.emotion.careersat)
getPar(RSA.emotion.careersat)
# plot RSA graph 
plot(RSA.emotion.careersat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1", "PA2"), 
     project = c("PA1", "PA2", "LOC"),        
     xlim    = c(-36,36),   
     ylim    = c(-36,36),    
     xlab = "Emotional Ability", 
     ylab = "Self View", 
     zlab = "Career Satisfaction"
)
# plot differently to see it better
RSA.emotion.careersat.reverse<-RSA(career.satisfaction~r.c.ERAtotal*r.c.ertselfrating,data=selfinsight.emo.careersat,
                                 center=FALSE,
                                 scale=FALSE,
                                 na.rm=TRUE,
                                 out.rm=NULL,
                                 breakline=FALSE, 
                                 models="default",
                                 cubic=FALSE,
                                 verbose=TRUE,
                                 estimator="MLR",
                                 se="robust",
                                 missing='listwise'
)
# plot reverse RSA graph 
plot(RSA.emotion.careersat.reverse,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-36,36),   
     ylim    = c(-36,36),    
     xlab = "Emotional Ability (reversed)", 
     ylab = "Self View (reversed)", 
     zlab = "career Satisfaction"
)

#### ROBUSTNESS CHECK WITH COMMON GRAND MEAN CENTERING
#### common grand mean centering
commonmean <- (mean(selfinsight.emo.careersat$ERAtotal)+mean(selfinsight.emo.careersat$ertselfrating))/2
selfinsight.emo.careersat$cc.ERAtotal <- (selfinsight.emo.careersat$ERAtotal)-commonmean
selfinsight.emo.careersat$cc.ertselfrating <- (selfinsight.emo.careersat$ertselfrating)-commonmean
with(selfinsight.emo.careersat,describe(cbind(cc.ERAtotal,ERAtotal,cc.ertselfrating,ertselfrating,commonmean)))
#### test for self-insight with polynomial regression analyses and RSA
RSA.emotion.careersat.cc<-RSA(career.satisfaction~cc.ERAtotal*cc.ertselfrating,data=selfinsight.emo.careersat,
                            center=FALSE,
                            scale=FALSE,
                            na.rm=TRUE,
                            out.rm=NULL,
                            breakline=FALSE, 
                            models="default",
                            cubic=FALSE,
                            verbose=TRUE,
                            estimator="MLR",
                            se="robust",
                            missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.emotion.careersat.cc)
getPar(RSA.emotion.careersat.cc)
#### plot RSA graph 
plot(RSA.emotion.careersat.cc,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-52.47,19.53),   
     ylim    = c(-52.47,19.53),    
     xlab = "Emotional Ability", 
     ylab = "Self View", 
     zlab = "Career Satisfaction"
)

#### ERA AND RELATIONSHIP SATISFACTION ####

#### calculate relationship satisfaction score for each diary

selfinsight$interpersonal_interpersonal_1_1 <- as.numeric(selfinsight$interpersonal_interpersonal_1_1); selfinsight$interpersonal_interpersonal_2_1 <- as.numeric(selfinsight$interpersonal_interpersonal_2_1); selfinsight$interpersonal_interpersonal_3_1 <- as.numeric(selfinsight$interpersonal_interpersonal_3_1); selfinsight$interpersonal_interpersonal_4_1 <- as.numeric(selfinsight$interpersonal_interpersonal_4_1)
selfinsight$interpersonal.diary1 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_1,interpersonal_interpersonal_2_1,interpersonal_interpersonal_3_1,interpersonal_interpersonal_4_1),na.rm=TRUE)))
selfinsight$interpersonal_interpersonal_1_2 <- as.numeric(selfinsight$interpersonal_interpersonal_1_2); selfinsight$interpersonal_interpersonal_2_2 <- as.numeric(selfinsight$interpersonal_interpersonal_2_2); selfinsight$interpersonal_interpersonal_3_2 <- as.numeric(selfinsight$interpersonal_interpersonal_3_2); selfinsight$interpersonal_interpersonal_4_2 <- as.numeric(selfinsight$interpersonal_interpersonal_4_2)
selfinsight$interpersonal.diary2 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_2,interpersonal_interpersonal_2_2,interpersonal_interpersonal_3_2,interpersonal_interpersonal_4_2),na.rm=TRUE)))
selfinsight$interpersonal_interpersonal_1_3 <- as.numeric(selfinsight$interpersonal_interpersonal_1_3); selfinsight$interpersonal_interpersonal_2_3 <- as.numeric(selfinsight$interpersonal_interpersonal_2_3); selfinsight$interpersonal_interpersonal_3_3 <- as.numeric(selfinsight$interpersonal_interpersonal_3_3); selfinsight$interpersonal_interpersonal_4_3 <- as.numeric(selfinsight$interpersonal_interpersonal_4_3)
selfinsight$interpersonal.diary3 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_3,interpersonal_interpersonal_2_3,interpersonal_interpersonal_3_3,interpersonal_interpersonal_4_3),na.rm=TRUE)))
selfinsight$interpersonal_interpersonal_1_4 <- as.numeric(selfinsight$interpersonal_interpersonal_1_4); selfinsight$interpersonal_interpersonal_2_4 <- as.numeric(selfinsight$interpersonal_interpersonal_2_4); selfinsight$interpersonal_interpersonal_3_4 <- as.numeric(selfinsight$interpersonal_interpersonal_3_4); selfinsight$interpersonal_interpersonal_4_4 <- as.numeric(selfinsight$interpersonal_interpersonal_4_4)
selfinsight$interpersonal.diary4 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_4,interpersonal_interpersonal_2_4,interpersonal_interpersonal_3_4,interpersonal_interpersonal_4_4),na.rm=TRUE)))
selfinsight$interpersonal_interpersonal_1_5 <- as.numeric(selfinsight$interpersonal_interpersonal_1_5); selfinsight$interpersonal_interpersonal_2_5 <- as.numeric(selfinsight$interpersonal_interpersonal_2_5); selfinsight$interpersonal_interpersonal_3_5 <- as.numeric(selfinsight$interpersonal_interpersonal_3_5); selfinsight$interpersonal_interpersonal_4_5 <- as.numeric(selfinsight$interpersonal_interpersonal_4_5)
selfinsight$interpersonal.diary5 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_5,interpersonal_interpersonal_2_5,interpersonal_interpersonal_3_5,interpersonal_interpersonal_4_5),na.rm=TRUE)))
selfinsight$interpersonal_interpersonal_1_6 <- as.numeric(selfinsight$interpersonal_interpersonal_1_6); selfinsight$interpersonal_interpersonal_2_6 <- as.numeric(selfinsight$interpersonal_interpersonal_2_6); selfinsight$interpersonal_interpersonal_3_6 <- as.numeric(selfinsight$interpersonal_interpersonal_3_6); selfinsight$interpersonal_interpersonal_4_6 <- as.numeric(selfinsight$interpersonal_interpersonal_4_6)
selfinsight$interpersonal.diary6 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_6,interpersonal_interpersonal_2_6,interpersonal_interpersonal_3_6,interpersonal_interpersonal_4_6),na.rm=TRUE)))
selfinsight$interpersonal_interpersonal_1_7 <- as.numeric(selfinsight$interpersonal_interpersonal_1_7); selfinsight$interpersonal_interpersonal_2_7 <- as.numeric(selfinsight$interpersonal_interpersonal_2_7); selfinsight$interpersonal_interpersonal_3_7 <- as.numeric(selfinsight$interpersonal_interpersonal_3_7); selfinsight$interpersonal_interpersonal_4_7 <- as.numeric(selfinsight$interpersonal_interpersonal_4_7)
selfinsight$interpersonal.diary7 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_7,interpersonal_interpersonal_2_7,interpersonal_interpersonal_3_7,interpersonal_interpersonal_4_7),na.rm=TRUE)))
selfinsight$interpersonal_interpersonal_1_8 <- as.numeric(selfinsight$interpersonal_interpersonal_1_8); selfinsight$interpersonal_interpersonal_2_8 <- as.numeric(selfinsight$interpersonal_interpersonal_2_8); selfinsight$interpersonal_interpersonal_3_8 <- as.numeric(selfinsight$interpersonal_interpersonal_3_8); selfinsight$interpersonal_interpersonal_4_8 <- as.numeric(selfinsight$interpersonal_interpersonal_4_8)
selfinsight$interpersonal.diary8 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_8,interpersonal_interpersonal_2_8,interpersonal_interpersonal_3_8,interpersonal_interpersonal_4_8),na.rm=TRUE)))
selfinsight$interpersonal_interpersonal_1_9 <- as.numeric(selfinsight$interpersonal_interpersonal_1_9); selfinsight$interpersonal_interpersonal_2_9 <- as.numeric(selfinsight$interpersonal_interpersonal_2_9); selfinsight$interpersonal_interpersonal_3_9 <- as.numeric(selfinsight$interpersonal_interpersonal_3_9); selfinsight$interpersonal_interpersonal_4_9 <- as.numeric(selfinsight$interpersonal_interpersonal_4_9)
selfinsight$interpersonal.diary9 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_9,interpersonal_interpersonal_2_9,interpersonal_interpersonal_3_9,interpersonal_interpersonal_4_9),na.rm=TRUE)))
selfinsight$interpersonal_interpersonal_1_0 <- as.numeric(selfinsight$interpersonal_interpersonal_1_0); selfinsight$interpersonal_interpersonal_2_0 <- as.numeric(selfinsight$interpersonal_interpersonal_2_0); selfinsight$interpersonal_interpersonal_3_0 <- as.numeric(selfinsight$interpersonal_interpersonal_3_0); selfinsight$interpersonal_interpersonal_4_0 <- as.numeric(selfinsight$interpersonal_interpersonal_4_0)
selfinsight$interpersonal.diary0 <- with(selfinsight,(rowMeans(cbind(interpersonal_interpersonal_1_0,interpersonal_interpersonal_2_0,interpersonal_interpersonal_3_0,interpersonal_interpersonal_4_0),na.rm=TRUE)))

# average scores for all diaries
selfinsight$relationship.satisfaction <- with(selfinsight,(rowMeans(cbind(interpersonal.diary1,interpersonal.diary2,interpersonal.diary3,interpersonal.diary4,interpersonal.diary5,interpersonal.diary6,interpersonal.diary7,interpersonal.diary8,interpersonal.diary9,interpersonal.diary0),na.rm=TRUE)))
with(selfinsight, describe(cbind(relationship.satisfaction,interpersonal.diary1,interpersonal.diary2,interpersonal.diary3,interpersonal.diary4,interpersonal.diary5,interpersonal.diary6,interpersonal.diary7,interpersonal.diary8,interpersonal.diary9,interpersonal.diary0)))

### REMOVE OUTLIERS (p. 32 of pre-registration)

### UNIVARIATE OUTLIERS
### observations 3.29 standard deviations above and below the mean (Tabachnick & Fidell)

meanERAtotal <- mean(selfinsight$ERAtotal); sdERAtotal <- SD(selfinsight$ERAtotal); lowcutoffERAtotal <- meanERAtotal-(3.29*sdERAtotal); highcutoffERAtotal <- meanERAtotal+(3.29*sdERAtotal)
table(selfinsight$ERAtotal)
meanERAtotal; sdERAtotal; lowcutoffERAtotal; highcutoffERAtotal 
#### note: 11 observations are LOW outliers on the emotion recognition test
#### as per the pre-registration, these analyses will be removed from the analyses
meanertselfrating <- mean(selfinsight$ertselfrating); sdertselfrating <- SD(selfinsight$ertselfrating); lowcutoffertselfrating <- meanertselfrating-(3.29*sdertselfrating); highcutoffertselfrating <- meanertselfrating+(3.29*sdertselfrating)
table(selfinsight$ertselfrating)
meanertselfrating; sdertselfrating; lowcutoffertselfrating; highcutoffertselfrating 
#### note: 15 observations are LOW outliers on self-views about emotion recognition
#### as per the pre-registration, these analyses will be removed from the analyses
meanrelationship.satisfaction <- mean(selfinsight$relationship.satisfaction); sdrelationship.satisfaction <- SD(selfinsight$relationship.satisfaction); lowcutoffrelationship.satisfaction <- meanrelationship.satisfaction-(3.29*sdrelationship.satisfaction); highcutoffrelationship.satisfaction <- meanrelationship.satisfaction+(3.29*sdrelationship.satisfaction)
table(selfinsight$relationship.satisfaction)
meanrelationship.satisfaction; sdrelationship.satisfaction; lowcutoffrelationship.satisfaction; highcutoffrelationship.satisfaction 
#### note: 2 observations are HIGH outliers on relationship satisfaction (LOW cut off is outside the possible range for the measure)
selfinsight.emo.relationshipsat.temp1 <- subset(selfinsight, ERAtotal<highcutoffERAtotal)
selfinsight.emo.relationshipsat.temp2 <- subset(selfinsight.emo.relationshipsat.temp1, ERAtotal>lowcutoffERAtotal)
selfinsight.emo.relationshipsat.temp3 <- subset(selfinsight.emo.relationshipsat.temp2, ertselfrating<highcutoffertselfrating)
selfinsight.emo.relationshipsat.temp4 <- subset(selfinsight.emo.relationshipsat.temp3, ertselfrating>lowcutoffertselfrating)
selfinsight.emo.relationshipsat.temp5 <- subset(selfinsight.emo.relationshipsat.temp4, relationship.satisfaction<highcutoffrelationship.satisfaction)
selfinsight.emo.relationshipsat <- subset(selfinsight.emo.relationshipsat.temp5, relationship.satisfaction>lowcutoffrelationship.satisfaction)

### multivariate outliers
### observations with values that exceed the cut-offs recommended by Bollen and Jackman98 and Belsley, Kuh, and Welsch on three indices: Cook’s D, Hat values, and difference in fits (DFFITS)
## leverage measures
selfinsight.emo.relationshipsat$out <- FALSE
selfinsight.emo.relationshipsat$c.ERAtotal2 <- (selfinsight.emo.relationshipsat$c.ERAtotal)^2
selfinsight.emo.relationshipsat$c.ertselfrating2 <- (selfinsight.emo.relationshipsat$c.ertselfrating)^2
selfinsight.emo.relationshipsat$c.ERAtotal_c.ertselfrating <- selfinsight.emo.relationshipsat$c.ERAtotal*selfinsight.emo.relationshipsat$c.ertselfrating
full_model <- lm(relationship.satisfaction~c.ERAtotal+c.ertselfrating+c.ERAtotal2+c.ERAtotal_c.ertselfrating+c.ertselfrating2, data=selfinsight.emo.relationshipsat, na.action=na.exclude); summary(full_model)
infl <- influence.measures(full_model)
# function to get leverage measures (Bollen & Jackman, 1985)
selfinsight.emo.relationshipsat$out <- apply(infl$is.inf[, c("dffit", "cook.d", "hat")], 1, sum) == 3
# if a case has significant cooks d, dffit, AND hat, then it gets flagged as an outlier
n.out <- sum(na.omit(selfinsight.emo.relationshipsat$out) == TRUE)
# count the number of outliers identified by the previous line
n.out
# no outliers in data
# if outliers existed, they would be taken out (via data=accuracy[selfinsight.emo.relationshipsat.temp6$out==FALSE, ])

#### multicollinearity
ols_vif_tol(full_model)

#### testing assumptions of OLS regression ####

# normality of residuals
ols_plot_resid_qq(full_model)
ols_test_normality(full_model)
ols_plot_resid_hist(full_model)

# heteroskedasticity
bptest(full_model)

#### the data violate assumptions of homoscedasticity; use robust standard errors

#### test of significance of polynomial terms
summary(full_model.3.a <- lm(relationship.satisfaction~c.ERAtotal+c.ertselfrating, data=selfinsight.emo.relationshipsat, na.action=na.exclude))
summary(full_model.3.b <- lm(relationship.satisfaction~c.ERAtotal+c.ertselfrating+c.ERAtotal2+c.ERAtotal_c.ertselfrating+c.ertselfrating2, data=selfinsight.emo.relationshipsat, na.action=na.exclude))
anova(full_model.3.a, full_model.3.b)

#### test for self-insight with polynomial regression analyses and RSA
RSA.emotion.relationshipsat<-RSA(relationship.satisfaction~c.ERAtotal*c.ertselfrating,data=selfinsight.emo.relationshipsat,
                           center=FALSE,
                           scale=FALSE,
                           na.rm=TRUE,
                           out.rm=NULL,
                           breakline=FALSE, 
                           models="default",
                           cubic=FALSE,
                           verbose=TRUE,
                           estimator="MLR",
                           se="robust",
                           missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.emotion.relationshipsat)
getPar(RSA.emotion.relationshipsat)
# plot RSA graph 
plot(RSA.emotion.relationshipsat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-36,36),   
     ylim    = c(-36,36),    
     xlab = "Emotional Ability", 
     ylab = "Self View", 
     zlab = "Relationship Satisfaction"
)

#### ROBUSTNESS CHECK WITH COMMON GRAND MEAN CENTERING

#### common grand mean centering
commonmean <- (mean(selfinsight.emo.relationshipsat$ERAtotal)+mean(selfinsight.emo.relationshipsat$ertselfrating))/2
selfinsight.emo.relationshipsat$cc.ERAtotal <- (selfinsight.emo.relationshipsat$ERAtotal)-commonmean
selfinsight.emo.relationshipsat$cc.ertselfrating <- (selfinsight.emo.relationshipsat$ertselfrating)-commonmean

with(selfinsight.emo.relationshipsat,describe(cbind(cc.ERAtotal,ERAtotal,cc.ertselfrating,ertselfrating,commonmean)))

#### test for self-insight with polynomial regression analyses and RSA
RSA.emotion.relationshipsat.cc<-RSA(relationship.satisfaction~cc.ERAtotal*cc.ertselfrating,data=selfinsight.emo.relationshipsat,
                              center=FALSE,
                              scale=FALSE,
                              na.rm=TRUE,
                              out.rm=NULL,
                              breakline=FALSE, 
                              models="default",
                              cubic=FALSE,
                              verbose=TRUE,
                              estimator="MLR",
                              se="robust",
                              missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.emotion.relationshipsat.cc)
getPar(RSA.emotion.relationshipsat.cc)
#### plot RSA graph 
plot(RSA.emotion.relationshipsat.cc,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-52.47,19.53),   
     ylim    = c(-52.47,19.53),    
     xlab = "Emotional Ability", 
     ylab = "Self View", 
     zlab = "Relationship Satisfaction"
)

#### ANALYSES WITH COGNITIVE ABILITY ####

#### CENTER VARIABLES

selfinsight.valid.cognitive.ability$c.raventotal <- selfinsight.valid.cognitive.ability$raventotal-7.5
selfinsight.valid.cognitive.ability$c.ravenselfrating <- selfinsight.valid.cognitive.ability$ravenselfrating-7.5
with(selfinsight.valid.cognitive.ability,describe(cbind(c.raventotal,raventotal,c.ravenselfrating,ravenselfrating)))

#### COGNITIVE ABILITY AND LIFE SATISFACTION ####

# calculate life satisfaction score for each diary in this data frame

selfinsight.valid.cognitive.ability$swls_swls_1_1 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_1); selfinsight.valid.cognitive.ability$swls_swls_2_1 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_1); selfinsight.valid.cognitive.ability$swls_swls_3_1 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_1); selfinsight.valid.cognitive.ability$swls_swls_4_1 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_1); selfinsight.valid.cognitive.ability$swls_swls_5_1 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_1)
selfinsight.valid.cognitive.ability$swls.diary1 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_1,swls_swls_2_1,swls_swls_3_1,swls_swls_4_1,swls_swls_5_1),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$swls_swls_1_2 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_2); selfinsight.valid.cognitive.ability$swls_swls_2_2 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_2); selfinsight.valid.cognitive.ability$swls_swls_3_2 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_2); selfinsight.valid.cognitive.ability$swls_swls_4_2 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_2); selfinsight.valid.cognitive.ability$swls_swls_5_2 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_2)
selfinsight.valid.cognitive.ability$swls.diary2 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_2,swls_swls_2_2,swls_swls_3_2,swls_swls_4_2,swls_swls_5_2),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$swls_swls_1_3 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_3); selfinsight.valid.cognitive.ability$swls_swls_2_3 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_3); selfinsight.valid.cognitive.ability$swls_swls_3_3 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_3); selfinsight.valid.cognitive.ability$swls_swls_4_3 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_3); selfinsight.valid.cognitive.ability$swls_swls_5_3 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_3)
selfinsight.valid.cognitive.ability$swls.diary3 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_3,swls_swls_2_3,swls_swls_3_3,swls_swls_4_3,swls_swls_5_3),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$swls_swls_1_4 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_4); selfinsight.valid.cognitive.ability$swls_swls_2_4 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_4); selfinsight.valid.cognitive.ability$swls_swls_3_4 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_4); selfinsight.valid.cognitive.ability$swls_swls_4_4 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_4); selfinsight.valid.cognitive.ability$swls_swls_5_4 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_4)
selfinsight.valid.cognitive.ability$swls.diary4 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_4,swls_swls_2_4,swls_swls_3_4,swls_swls_4_4,swls_swls_5_4),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$swls_swls_1_5 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_5); selfinsight.valid.cognitive.ability$swls_swls_2_5 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_5); selfinsight.valid.cognitive.ability$swls_swls_3_5 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_5); selfinsight.valid.cognitive.ability$swls_swls_4_5 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_5); selfinsight.valid.cognitive.ability$swls_swls_5_5 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_5)
selfinsight.valid.cognitive.ability$swls.diary5 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_5,swls_swls_2_5,swls_swls_3_5,swls_swls_4_5,swls_swls_5_5),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$swls_swls_1_6 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_6); selfinsight.valid.cognitive.ability$swls_swls_2_6 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_6); selfinsight.valid.cognitive.ability$swls_swls_3_6 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_6); selfinsight.valid.cognitive.ability$swls_swls_4_6 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_6); selfinsight.valid.cognitive.ability$swls_swls_5_6 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_6)
selfinsight.valid.cognitive.ability$swls.diary6 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_6,swls_swls_2_6,swls_swls_3_6,swls_swls_4_6,swls_swls_5_6),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$swls_swls_1_7 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_7); selfinsight.valid.cognitive.ability$swls_swls_2_7 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_7); selfinsight.valid.cognitive.ability$swls_swls_3_7 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_7); selfinsight.valid.cognitive.ability$swls_swls_4_7 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_7); selfinsight.valid.cognitive.ability$swls_swls_5_7 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_7)
selfinsight.valid.cognitive.ability$swls.diary7 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_7,swls_swls_2_7,swls_swls_3_7,swls_swls_4_7,swls_swls_5_7),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$swls_swls_1_8 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_8); selfinsight.valid.cognitive.ability$swls_swls_2_8 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_8); selfinsight.valid.cognitive.ability$swls_swls_3_8 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_8); selfinsight.valid.cognitive.ability$swls_swls_4_8 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_8); selfinsight.valid.cognitive.ability$swls_swls_5_8 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_8)
selfinsight.valid.cognitive.ability$swls.diary8 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_8,swls_swls_2_8,swls_swls_3_8,swls_swls_4_8,swls_swls_5_8),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$swls_swls_1_9 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_9); selfinsight.valid.cognitive.ability$swls_swls_2_9 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_9); selfinsight.valid.cognitive.ability$swls_swls_3_9 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_9); selfinsight.valid.cognitive.ability$swls_swls_4_9 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_9); selfinsight.valid.cognitive.ability$swls_swls_5_9 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_9)
selfinsight.valid.cognitive.ability$swls.diary9 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_9,swls_swls_2_9,swls_swls_3_9,swls_swls_4_9,swls_swls_5_9),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$swls_swls_1_0 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_1_0); selfinsight.valid.cognitive.ability$swls_swls_2_0 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_2_0); selfinsight.valid.cognitive.ability$swls_swls_3_0 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_3_0); selfinsight.valid.cognitive.ability$swls_swls_4_0 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_4_0); selfinsight.valid.cognitive.ability$swls_swls_5_0 <- as.numeric(selfinsight.valid.cognitive.ability$swls_swls_5_0)
selfinsight.valid.cognitive.ability$swls.diary0 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls_swls_1_0,swls_swls_2_0,swls_swls_3_0,swls_swls_4_0,swls_swls_5_0),na.rm=TRUE)))

# average scores for all diaries
selfinsight.valid.cognitive.ability$life.satisfaction <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(swls.diary1,swls.diary2,swls.diary3,swls.diary4,swls.diary5,swls.diary6,swls.diary7,swls.diary8,swls.diary9,swls.diary0),na.rm=TRUE)))
with(selfinsight.valid.cognitive.ability, describe(cbind(life.satisfaction,swls.diary1,swls.diary2,swls.diary3,swls.diary4,swls.diary5,swls.diary6,swls.diary7,swls.diary8,swls.diary9,swls.diary0)))

#### REMOVE OUTLIERS (p. 32 of pre-registration)

#### UNIVARIATE OUTLIERS
#### observations 3.29 standard deviations above and below the mean (Tabachnick & Fidell)
meanraventotal <- mean(selfinsight.valid.cognitive.ability$raventotal); sdraventotal <- SD(selfinsight.valid.cognitive.ability$raventotal); lowcutoffraventotal <- meanraventotal-(3.29*sdraventotal); highcutoffraventotal <- meanraventotal+(3.29*sdraventotal)
table(selfinsight.valid.cognitive.ability$raventotal)
meanraventotal; sdraventotal; lowcutoffraventotal; highcutoffraventotal 
#### note: NO observations are outliers on the cognitive ability test (cut offs are outside the range)
meanravenselfrating <- mean(selfinsight.valid.cognitive.ability$ravenselfrating); sdravenselfrating <- SD(selfinsight.valid.cognitive.ability$ravenselfrating); lowcutoffravenselfrating <- meanravenselfrating-(3.29*sdravenselfrating); highcutoffravenselfrating <- meanravenselfrating+(3.29*sdravenselfrating)
table(selfinsight.valid.cognitive.ability$ravenselfrating)
meanravenselfrating; sdravenselfrating; lowcutoffravenselfrating; highcutoffravenselfrating 
#### note: NO observations are outliers on self-views about cognitive ability (cut offs are outside the range)
meanlife.satisfaction <- mean(selfinsight.valid.cognitive.ability$life.satisfaction); sdlife.satisfaction <- SD(selfinsight.valid.cognitive.ability$life.satisfaction); lowcutofflife.satisfaction <- meanlife.satisfaction-(3.29*sdlife.satisfaction); highcutofflife.satisfaction <- meanlife.satisfaction+(3.29*sdlife.satisfaction)
meanlife.satisfaction; sdlife.satisfaction; lowcutofflife.satisfaction; highcutofflife.satisfaction 
#### note: NO observations are outliers on life satisfaction (cut offs are actually outside the possible range for the measure)
selfinsight.cog.lifesat.temp1 <- subset(selfinsight.valid.cognitive.ability, raventotal<highcutoffraventotal)
selfinsight.cog.lifesat.temp2 <- subset(selfinsight.cog.lifesat.temp1, raventotal>lowcutoffraventotal)
selfinsight.cog.lifesat.temp3 <- subset(selfinsight.cog.lifesat.temp2, ravenselfrating<highcutoffravenselfrating)
selfinsight.cog.lifesat.temp4 <- subset(selfinsight.cog.lifesat.temp3, ravenselfrating>lowcutoffravenselfrating)
selfinsight.cog.lifesat.temp5 <- subset(selfinsight.cog.lifesat.temp4, life.satisfaction<highcutofflife.satisfaction)
selfinsight.cog.lifesat <- subset(selfinsight.cog.lifesat.temp5, life.satisfaction>lowcutofflife.satisfaction)

#### multivariate outliers
#### observations with values that exceed the cut-offs recommended by Bollen and Jackman98 and Belsley, Kuh, and Welsch on three indices: Cook’s D, Hat values, and difference in fits (DFFITS)
#### leverage measures
selfinsight.cog.lifesat$out <- FALSE
selfinsight.cog.lifesat$c.raventotal2 <- (selfinsight.cog.lifesat$c.raventotal)^2
selfinsight.cog.lifesat$c.ravenselfrating2 <- (selfinsight.cog.lifesat$c.ravenselfrating)^2
selfinsight.cog.lifesat$c.raventotal_c.ravenselfrating <- selfinsight.cog.lifesat$c.raventotal*selfinsight.cog.lifesat$c.ravenselfrating
full_model <- lm(life.satisfaction~c.raventotal+c.ravenselfrating+c.raventotal2+c.raventotal_c.ravenselfrating+c.ravenselfrating2, data=selfinsight.cog.lifesat, na.action=na.exclude); summary(full_model)
infl <- influence.measures(full_model)
# function to get leverage measures (Bollen & Jackman, 1985)
selfinsight.cog.lifesat$out <- apply(infl$is.inf[, c("dffit", "cook.d", "hat")], 1, sum) == 3
# if a case has significant cooks d, dffit, AND hat, then it gets flagged as an outlier
n.out <- sum(na.omit(selfinsight.cog.lifesat$out) == TRUE)
# count the number of outliers identified by the previous line
n.out
# no outliers in data
# if outliers existed, they would be taken out (via data=accuracy[selfinsight.cog.lifesat.temp6$out==FALSE, ])

#### multicollinearity
ols_vif_tol(full_model)

#### testing assumptions of OLS regression ####

# normality of residuals
ols_plot_resid_qq(full_model)
ols_test_normality(full_model)
ols_plot_resid_hist(full_model)

# heteroskedasticity
bptest(full_model)

#### the data violate assumptions of homoscedasticity; use robust standard errors

#### NO test of significance of polynomial terms because full model is not significant

#### test for self-insight with polynomial regression analyses and RSA
RSA.cognition.lifesat<-RSA(life.satisfaction~c.raventotal*c.ravenselfrating,data=selfinsight.cog.lifesat,
                         center=FALSE,
                         scale=FALSE,
                         na.rm=TRUE,
                         out.rm=NULL,
                         breakline=FALSE, 
                         models="default",
                         cubic=FALSE,
                         verbose=TRUE,
                         estimator="MLR",
                         se="robust",
                         missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.cognition.lifesat)
getPar(RSA.cognition.lifesat)
# plot RSA graph 
plot(RSA.cognition.lifesat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-7.5,7.5),   
     ylim    = c(-7.5,7.5),
     zlim    = c(1,7),
     xlab = "Cognitive Ability", 
     ylab = "Self View", 
     zlab = "Life Satisfaction"
)

#### ROBUSTNESS CHECK WITH COMMON GRAND MEAN CENTERING

#### common grand mean centering
commonmean <- (mean(selfinsight.cog.lifesat$raventotal)+mean(selfinsight.cog.lifesat$ravenselfrating))/2
selfinsight.cog.lifesat$cc.raventotal <- (selfinsight.cog.lifesat$raventotal)-commonmean
selfinsight.cog.lifesat$cc.ravenselfrating <- (selfinsight.cog.lifesat$ravenselfrating)-commonmean

with(selfinsight.cog.lifesat,describe(cbind(cc.raventotal,raventotal,cc.ravenselfrating,ravenselfrating,commonmean)))

#### test for self-insight with polynomial regression analyses and RSA
RSA.cognition.lifesat.cc<-RSA(life.satisfaction~cc.raventotal*cc.ravenselfrating,data=selfinsight.cog.lifesat,
                            center=FALSE,
                            scale=FALSE,
                            na.rm=TRUE,
                            out.rm=NULL,
                            breakline=FALSE, 
                            models="default",
                            cubic=FALSE,
                            verbose=TRUE,
                            estimator="MLR",
                            se="robust",
                            missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.cognition.lifesat.cc)
getPar(RSA.cognition.lifesat.cc)
# plot RSA graph 
plot(RSA.cognition.lifesat.cc,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-6.95,8.05),   
     ylim    = c(-6.95,8.05),    
     zlim    = c(1,7),
     xlab = "Cognitive Ability", 
     ylab = "Self View", 
     zlab = "Life Satisfaction"
)

#### COGNITIVE ABILITY AND CAREER SATISFACTION ####

#### calculate career satisfaction score for each diary in this data frame

selfinsight.valid.cognitive.ability$institutional_institutional_1_1 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_1); selfinsight.valid.cognitive.ability$institutional_institutional_2_1 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_1); selfinsight.valid.cognitive.ability$institutional_institutional_3_1 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_1); selfinsight.valid.cognitive.ability$institutional_institutional_4_1 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_1); selfinsight.valid.cognitive.ability$institutional_institutional_5_1 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_1)
selfinsight.valid.cognitive.ability$institutional.diary1 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_1,institutional_institutional_2_1,institutional_institutional_3_1,institutional_institutional_4_1,institutional_institutional_5_1),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$institutional_institutional_1_2 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_2); selfinsight.valid.cognitive.ability$institutional_institutional_2_2 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_2); selfinsight.valid.cognitive.ability$institutional_institutional_3_2 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_2); selfinsight.valid.cognitive.ability$institutional_institutional_4_2 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_2); selfinsight.valid.cognitive.ability$institutional_institutional_5_2 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_2)
selfinsight.valid.cognitive.ability$institutional.diary2 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_2,institutional_institutional_2_2,institutional_institutional_3_2,institutional_institutional_4_2,institutional_institutional_5_2),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$institutional_institutional_1_3 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_3); selfinsight.valid.cognitive.ability$institutional_institutional_2_3 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_3); selfinsight.valid.cognitive.ability$institutional_institutional_3_3 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_3); selfinsight.valid.cognitive.ability$institutional_institutional_4_3 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_3); selfinsight.valid.cognitive.ability$institutional_institutional_5_3 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_3)
selfinsight.valid.cognitive.ability$institutional.diary3 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_3,institutional_institutional_2_3,institutional_institutional_3_3,institutional_institutional_4_3,institutional_institutional_5_3),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$institutional_institutional_1_4 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_4); selfinsight.valid.cognitive.ability$institutional_institutional_2_4 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_4); selfinsight.valid.cognitive.ability$institutional_institutional_3_4 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_4); selfinsight.valid.cognitive.ability$institutional_institutional_4_4 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_4); selfinsight.valid.cognitive.ability$institutional_institutional_5_4 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_4)
selfinsight.valid.cognitive.ability$institutional.diary4 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_4,institutional_institutional_2_4,institutional_institutional_3_4,institutional_institutional_4_4,institutional_institutional_5_4),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$institutional_institutional_1_5 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_5); selfinsight.valid.cognitive.ability$institutional_institutional_2_5 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_5); selfinsight.valid.cognitive.ability$institutional_institutional_3_5 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_5); selfinsight.valid.cognitive.ability$institutional_institutional_4_5 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_5); selfinsight.valid.cognitive.ability$institutional_institutional_5_5 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_5)
selfinsight.valid.cognitive.ability$institutional.diary5 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_5,institutional_institutional_2_5,institutional_institutional_3_5,institutional_institutional_4_5,institutional_institutional_5_5),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$institutional_institutional_1_6 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_6); selfinsight.valid.cognitive.ability$institutional_institutional_2_6 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_6); selfinsight.valid.cognitive.ability$institutional_institutional_3_6 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_6); selfinsight.valid.cognitive.ability$institutional_institutional_4_6 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_6); selfinsight.valid.cognitive.ability$institutional_institutional_5_6 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_6)
selfinsight.valid.cognitive.ability$institutional.diary6 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_6,institutional_institutional_2_6,institutional_institutional_3_6,institutional_institutional_4_6,institutional_institutional_5_6),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$institutional_institutional_1_7 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_7); selfinsight.valid.cognitive.ability$institutional_institutional_2_7 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_7); selfinsight.valid.cognitive.ability$institutional_institutional_3_7 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_7); selfinsight.valid.cognitive.ability$institutional_institutional_4_7 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_7); selfinsight.valid.cognitive.ability$institutional_institutional_5_7 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_7)
selfinsight.valid.cognitive.ability$institutional.diary7 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_7,institutional_institutional_2_7,institutional_institutional_3_7,institutional_institutional_4_7,institutional_institutional_5_7),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$institutional_institutional_1_8 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_8); selfinsight.valid.cognitive.ability$institutional_institutional_2_8 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_8); selfinsight.valid.cognitive.ability$institutional_institutional_3_8 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_8); selfinsight.valid.cognitive.ability$institutional_institutional_4_8 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_8); selfinsight.valid.cognitive.ability$institutional_institutional_5_8 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_8)
selfinsight.valid.cognitive.ability$institutional.diary8 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_8,institutional_institutional_2_8,institutional_institutional_3_8,institutional_institutional_4_8,institutional_institutional_5_8),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$institutional_institutional_1_9 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_9); selfinsight.valid.cognitive.ability$institutional_institutional_2_9 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_9); selfinsight.valid.cognitive.ability$institutional_institutional_3_9 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_9); selfinsight.valid.cognitive.ability$institutional_institutional_4_9 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_9); selfinsight.valid.cognitive.ability$institutional_institutional_5_9 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_9)
selfinsight.valid.cognitive.ability$institutional.diary9 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_9,institutional_institutional_2_9,institutional_institutional_3_9,institutional_institutional_4_9,institutional_institutional_5_9),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$institutional_institutional_1_0 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_1_0); selfinsight.valid.cognitive.ability$institutional_institutional_2_0 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_2_0); selfinsight.valid.cognitive.ability$institutional_institutional_3_0 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_3_0); selfinsight.valid.cognitive.ability$institutional_institutional_4_0 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_4_0); selfinsight.valid.cognitive.ability$institutional_institutional_5_0 <- as.numeric(selfinsight.valid.cognitive.ability$institutional_institutional_5_0)
selfinsight.valid.cognitive.ability$institutional.diary0 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional_institutional_1_0,institutional_institutional_2_0,institutional_institutional_3_0,institutional_institutional_4_0,institutional_institutional_5_0),na.rm=TRUE)))

# average scores for all diaries
selfinsight.valid.cognitive.ability$career.satisfaction <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(institutional.diary1,institutional.diary2,institutional.diary3,institutional.diary4,institutional.diary5,institutional.diary6,institutional.diary7,institutional.diary8,institutional.diary9,institutional.diary0),na.rm=TRUE)))
with(selfinsight.valid.cognitive.ability, describe(cbind(career.satisfaction,institutional.diary1,institutional.diary2,institutional.diary3,institutional.diary4,institutional.diary5,institutional.diary6,institutional.diary7,institutional.diary8,institutional.diary9,institutional.diary0)))

#### REMOVE OUTLIERS (p. 32 of pre-registration)

#### UNIVARIATE OUTLIERS
#### observations 3.29 standard deviations above and below the mean (Tabachnick & Fidell)

meanraventotal <- mean(selfinsight.valid.cognitive.ability$raventotal); sdraventotal <- SD(selfinsight.valid.cognitive.ability$raventotal); lowcutoffraventotal <- meanraventotal-(3.29*sdraventotal); highcutoffraventotal <- meanraventotal+(3.29*sdraventotal)
table(selfinsight.valid.cognitive.ability$raventotal)
meanraventotal; sdraventotal; lowcutoffraventotal; highcutoffraventotal 
#### note: NO observations are outliers on the cognitive ability test (cut offs are outside the range)
meanravenselfrating <- mean(selfinsight.valid.cognitive.ability$ravenselfrating); sdravenselfrating <- SD(selfinsight.valid.cognitive.ability$ravenselfrating); lowcutoffravenselfrating <- meanravenselfrating-(3.29*sdravenselfrating); highcutoffravenselfrating <- meanravenselfrating+(3.29*sdravenselfrating)
table(selfinsight.valid.cognitive.ability$ravenselfrating)
meanravenselfrating; sdravenselfrating; lowcutoffravenselfrating; highcutoffravenselfrating 
#### note: NO observations are outliers on self-views about cognitive ability (cut offs are outside the range)
meancareer.satisfaction <- mean(selfinsight.valid.cognitive.ability$career.satisfaction); sdcareer.satisfaction <- SD(selfinsight.valid.cognitive.ability$career.satisfaction); lowcutoffcareer.satisfaction <- meancareer.satisfaction-(3.29*sdcareer.satisfaction); highcutoffcareer.satisfaction <- meancareer.satisfaction+(3.29*sdcareer.satisfaction)
meancareer.satisfaction; sdcareer.satisfaction; lowcutoffcareer.satisfaction; highcutoffcareer.satisfaction 
#### note: NO observations are outliers on career satisfaction (cut offs are actually outside the possible range for the measure)
selfinsight.cog.careersat.temp1 <- subset(selfinsight.valid.cognitive.ability, raventotal<highcutoffraventotal)
selfinsight.cog.careersat.temp2 <- subset(selfinsight.cog.careersat.temp1, raventotal>lowcutoffraventotal)
selfinsight.cog.careersat.temp3 <- subset(selfinsight.cog.careersat.temp2, ravenselfrating<highcutoffravenselfrating)
selfinsight.cog.careersat.temp4 <- subset(selfinsight.cog.careersat.temp3, ravenselfrating>lowcutoffravenselfrating)
selfinsight.cog.careersat.temp5 <- subset(selfinsight.cog.careersat.temp4, career.satisfaction<highcutoffcareer.satisfaction)
selfinsight.cog.careersat <- subset(selfinsight.cog.careersat.temp5, career.satisfaction>lowcutoffcareer.satisfaction)

#### multivariate outliers
#### observations with values that exceed the cut-offs recommended by Bollen and Jackman98 and Belsley, Kuh, and Welsch on three indices: Cook’s D, Hat values, and difference in fits (DFFITS)
#### leverage measures
selfinsight.cog.careersat$out <- FALSE
selfinsight.cog.careersat$c.raventotal2 <- (selfinsight.cog.careersat$c.raventotal)^2
selfinsight.cog.careersat$c.ravenselfrating2 <- (selfinsight.cog.careersat$c.ravenselfrating)^2
selfinsight.cog.careersat$c.raventotal_c.ravenselfrating <- selfinsight.cog.careersat$c.raventotal*selfinsight.cog.careersat$c.ravenselfrating
full_model <- lm(career.satisfaction~c.raventotal+c.ravenselfrating+c.raventotal2+c.raventotal_c.ravenselfrating+c.ravenselfrating2, data=selfinsight.cog.careersat, na.action=na.exclude); summary(full_model)
infl <- influence.measures(full_model)
# function to get leverage measures (Bollen & Jackman, 1985)
selfinsight.cog.careersat$out <- apply(infl$is.inf[, c("dffit", "cook.d", "hat")], 1, sum) == 3
# if a case has significant cooks d, dffit, AND hat, then it gets flagged as an outlier
n.out <- sum(na.omit(selfinsight.cog.careersat$out) == TRUE)
# count the number of outliers identified by the previous line
n.out
# no outliers in data
# if outliers existed, they would be taken out (via data=accuracy[selfinsight.cog.careersat.temp6$out==FALSE, ])

#### multicollinearity
ols_vif_tol(full_model)

#### testing assumptions of OLS regression ####

# normality of residuals
ols_plot_resid_qq(full_model)
ols_test_normality(full_model)
ols_plot_resid_hist(full_model)

# heteroskedasticity
bptest(full_model)

#### the data do not violate assumptions of homoscedasticity

#### test of significance of polynomial terms
summary(full_model.5.a <- lm(career.satisfaction~c.raventotal+c.ravenselfrating, data=selfinsight.cog.careersat, na.action=na.exclude))
summary(full_model.5.b <- lm(career.satisfaction~c.raventotal+c.ravenselfrating+c.raventotal2+c.raventotal_c.ravenselfrating+c.ravenselfrating2, data=selfinsight.cog.careersat, na.action=na.exclude))
anova(full_model.5.a, full_model.5.b)

#### test for self-insight with polynomial regression analyses and RSA
RSA.cognition.careersat<-RSA(career.satisfaction~c.raventotal*c.ravenselfrating,data=selfinsight.cog.careersat,
                           center=FALSE,
                           scale=FALSE,
                           na.rm=TRUE,
                           out.rm=NULL,
                           breakline=FALSE, 
                           models="default",
                           cubic=FALSE,
                           verbose=TRUE,
                           estimator="MLR",
                           se="robust",
                           missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.cognition.careersat)
getPar(RSA.cognition.careersat)
# plot RSA graph 
plot(RSA.cognition.careersat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-7.5,7.5),   
     ylim    = c(-7.5,7.5),
     zlim    = c(1,5),
     xlab = "Cognitive Ability", 
     ylab = "Self View", 
     zlab = "Career Satisfaction"
)

#### ROBUSTNESS CHECK WITH COMMON GRAND MEAN CENTERING
#### common grand mean centering
commonmean <- (mean(selfinsight.cog.careersat$raventotal)+mean(selfinsight.cog.careersat$ravenselfrating))/2
selfinsight.cog.careersat$cc.raventotal <- (selfinsight.cog.careersat$raventotal)-commonmean
selfinsight.cog.careersat$cc.ravenselfrating <- (selfinsight.cog.careersat$ravenselfrating)-commonmean
with(selfinsight.cog.careersat,describe(cbind(cc.raventotal,raventotal,cc.ravenselfrating,ravenselfrating,commonmean)))
#### test for self-insight with polynomial regression analyses and RSA
RSA.cognition.careersat.cc<-RSA(career.satisfaction~cc.raventotal*cc.ravenselfrating,data=selfinsight.cog.careersat,
                              center=FALSE,
                              scale=FALSE,
                              na.rm=TRUE,
                              out.rm=NULL,
                              breakline=FALSE, 
                              models="default",
                              cubic=FALSE,
                              verbose=TRUE,
                              estimator="MLR",
                              se="robust",
                              missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.cognition.careersat.cc)
getPar(RSA.cognition.careersat.cc)
# plot RSA graph 
plot(RSA.cognition.careersat.cc,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-10,10),   
     ylim    = c(-10,10),    
     zlim    = c(1,5),
     xlab = "Cognitive Ability", 
     ylab = "Self View", 
     zlab = "Career Satisfaction"
)

#### COGNITIVE ABILITY AND RELATIONSHIP SATISFACTION ####

#### calculate relationship satisfaction score for each diary in this data frame

selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_1 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_1); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_1 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_1); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_1 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_1); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_1 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_1)
selfinsight.valid.cognitive.ability$interpersonal.diary1 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_1,interpersonal_interpersonal_2_1,interpersonal_interpersonal_3_1,interpersonal_interpersonal_4_1),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_2 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_2); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_2 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_2); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_2 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_2); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_2 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_2)
selfinsight.valid.cognitive.ability$interpersonal.diary2 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_2,interpersonal_interpersonal_2_2,interpersonal_interpersonal_3_2,interpersonal_interpersonal_4_2),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_3 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_3); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_3 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_3); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_3 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_3); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_3 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_3)
selfinsight.valid.cognitive.ability$interpersonal.diary3 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_3,interpersonal_interpersonal_2_3,interpersonal_interpersonal_3_3,interpersonal_interpersonal_4_3),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_4 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_4); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_4 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_4); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_4 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_4); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_4 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_4)
selfinsight.valid.cognitive.ability$interpersonal.diary4 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_4,interpersonal_interpersonal_2_4,interpersonal_interpersonal_3_4,interpersonal_interpersonal_4_4),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_5 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_5); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_5 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_5); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_5 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_5); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_5 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_5)
selfinsight.valid.cognitive.ability$interpersonal.diary5 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_5,interpersonal_interpersonal_2_5,interpersonal_interpersonal_3_5,interpersonal_interpersonal_4_5),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_6 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_6); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_6 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_6); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_6 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_6); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_6 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_6)
selfinsight.valid.cognitive.ability$interpersonal.diary6 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_6,interpersonal_interpersonal_2_6,interpersonal_interpersonal_3_6,interpersonal_interpersonal_4_6),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_7 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_7); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_7 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_7); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_7 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_7); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_7 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_7)
selfinsight.valid.cognitive.ability$interpersonal.diary7 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_7,interpersonal_interpersonal_2_7,interpersonal_interpersonal_3_7,interpersonal_interpersonal_4_7),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_8 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_8); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_8 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_8); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_8 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_8); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_8 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_8)
selfinsight.valid.cognitive.ability$interpersonal.diary8 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_8,interpersonal_interpersonal_2_8,interpersonal_interpersonal_3_8,interpersonal_interpersonal_4_8),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_9 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_9); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_9 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_9); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_9 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_9); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_9 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_9)
selfinsight.valid.cognitive.ability$interpersonal.diary9 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_9,interpersonal_interpersonal_2_9,interpersonal_interpersonal_3_9,interpersonal_interpersonal_4_9),na.rm=TRUE)))
selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_0 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_1_0); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_0 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_2_0); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_0 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_3_0); selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_0 <- as.numeric(selfinsight.valid.cognitive.ability$interpersonal_interpersonal_4_0)
selfinsight.valid.cognitive.ability$interpersonal.diary0 <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal_interpersonal_1_0,interpersonal_interpersonal_2_0,interpersonal_interpersonal_3_0,interpersonal_interpersonal_4_0),na.rm=TRUE)))

# average scores for all diaries
selfinsight.valid.cognitive.ability$relationship.satisfaction <- with(selfinsight.valid.cognitive.ability,(rowMeans(cbind(interpersonal.diary1,interpersonal.diary2,interpersonal.diary3,interpersonal.diary4,interpersonal.diary5,interpersonal.diary6,interpersonal.diary7,interpersonal.diary8,interpersonal.diary9,interpersonal.diary0),na.rm=TRUE)))
with(selfinsight.valid.cognitive.ability, describe(cbind(relationship.satisfaction,interpersonal.diary1,interpersonal.diary2,interpersonal.diary3,interpersonal.diary4,interpersonal.diary5,interpersonal.diary6,interpersonal.diary7,interpersonal.diary8,interpersonal.diary9,interpersonal.diary0)))

#### REMOVE OUTLIERS (p. 32 of pre-registration)

#### UNIVARIATE OUTLIERS
#### observations 3.29 standard deviations above and below the mean (Tabachnick & Fidell)

meanraventotal <- mean(selfinsight.valid.cognitive.ability$raventotal); sdraventotal <- SD(selfinsight.valid.cognitive.ability$raventotal); lowcutoffraventotal <- meanraventotal-(3.29*sdraventotal); highcutoffraventotal <- meanraventotal+(3.29*sdraventotal)
table(selfinsight.valid.cognitive.ability$raventotal)
meanraventotal; sdraventotal; lowcutoffraventotal; highcutoffraventotal 
#### note: NO observations are outliers on the cognitive ability test (cut offs are outside the range)
meanravenselfrating <- mean(selfinsight.valid.cognitive.ability$ravenselfrating); sdravenselfrating <- SD(selfinsight.valid.cognitive.ability$ravenselfrating); lowcutoffravenselfrating <- meanravenselfrating-(3.29*sdravenselfrating); highcutoffravenselfrating <- meanravenselfrating+(3.29*sdravenselfrating)
table(selfinsight.valid.cognitive.ability$ravenselfrating)
meanravenselfrating; sdravenselfrating; lowcutoffravenselfrating; highcutoffravenselfrating 
#### note: NO observations are outliers on self-views about cognitive ability (cut offs are outside the range)
meanrelationship.satisfaction <- mean(selfinsight.valid.cognitive.ability$relationship.satisfaction); sdrelationship.satisfaction <- SD(selfinsight.valid.cognitive.ability$relationship.satisfaction); lowcutoffrelationship.satisfaction <- meanrelationship.satisfaction-(3.29*sdrelationship.satisfaction); highcutoffrelationship.satisfaction <- meanrelationship.satisfaction+(3.29*sdrelationship.satisfaction)
table(selfinsight.valid.cognitive.ability$relationship.satisfaction)
meanrelationship.satisfaction; sdrelationship.satisfaction; lowcutoffrelationship.satisfaction; highcutoffrelationship.satisfaction 
#### note: 2 observations are HIGH outliers on relationship satisfaction (LOW cut off is outside the possible range for the measure)
selfinsight.cog.relationshipsat.temp1 <- subset(selfinsight.valid.cognitive.ability, raventotal<highcutoffraventotal)
selfinsight.cog.relationshipsat.temp2 <- subset(selfinsight.cog.relationshipsat.temp1, raventotal>lowcutoffraventotal)
selfinsight.cog.relationshipsat.temp3 <- subset(selfinsight.cog.relationshipsat.temp2, ravenselfrating<highcutoffravenselfrating)
selfinsight.cog.relationshipsat.temp4 <- subset(selfinsight.cog.relationshipsat.temp3, ravenselfrating>lowcutoffravenselfrating)
selfinsight.cog.relationshipsat.temp5 <- subset(selfinsight.cog.relationshipsat.temp4, relationship.satisfaction<highcutoffrelationship.satisfaction)
selfinsight.cog.relationshipsat <- subset(selfinsight.cog.relationshipsat.temp5, relationship.satisfaction>lowcutoffrelationship.satisfaction)

#### multivariate outliers
#### observations with values that exceed the cut-offs recommended by Bollen and Jackman98 and Belsley, Kuh, and Welsch on three indices: Cook’s D, Hat values, and difference in fits (DFFITS)
#### leverage measures
selfinsight.cog.relationshipsat$out <- FALSE
selfinsight.cog.relationshipsat$c.raventotal2 <- (selfinsight.cog.relationshipsat$c.raventotal)^2
selfinsight.cog.relationshipsat$c.ravenselfrating2 <- (selfinsight.cog.relationshipsat$c.ravenselfrating)^2
selfinsight.cog.relationshipsat$c.raventotal_c.ravenselfrating <- selfinsight.cog.relationshipsat$c.raventotal*selfinsight.cog.relationshipsat$c.ravenselfrating
full_model <- lm(relationship.satisfaction~c.raventotal+c.ravenselfrating+c.raventotal2+c.raventotal_c.ravenselfrating+c.ravenselfrating2, data=selfinsight.cog.relationshipsat, na.action=na.exclude); summary(full_model)
infl <- influence.measures(full_model)
# function to get leverage measures (Bollen & Jackman, 1985)
selfinsight.cog.relationshipsat$out <- apply(infl$is.inf[, c("dffit", "cook.d", "hat")], 1, sum) == 3
# if a case has significant cooks d, dffit, AND hat, then it gets flagged as an outlier
n.out <- sum(na.omit(selfinsight.cog.relationshipsat$out) == TRUE)
# count the number of outliers identified by the previous line
n.out
# no outliers in data
# if outliers existed, they would be taken out (via data=accuracy[selfinsight.cog.relationshipsat.temp6$out==FALSE, ])

#### multicollinearity
ols_vif_tol(full_model)

#### testing assumptions of OLS regression ####

# normality of residuals
ols_plot_resid_qq(full_model)
ols_test_normality(full_model)
ols_plot_resid_hist(full_model)

# heteroskedasticity
bptest(full_model)

#### the data do not violate assumptions of homoscedasticity

#### test of significance of polynomial terms
summary(full_model.6.a <- lm(relationship.satisfaction~c.raventotal+c.ravenselfrating, data=selfinsight.cog.relationshipsat, na.action=na.exclude))
summary(full_model.6.b <- lm(relationship.satisfaction~c.raventotal+c.ravenselfrating+c.raventotal2+c.raventotal_c.ravenselfrating+c.ravenselfrating2, data=selfinsight.cog.relationshipsat, na.action=na.exclude))
anova(full_model.6.a, full_model.6.b)


#### test for self-insight with polynomial regression analyses and RSA
RSA.cognition.relationshipsat<-RSA(relationship.satisfaction~c.raventotal*c.ravenselfrating,data=selfinsight.cog.relationshipsat,
                                 center=FALSE,
                                 scale=FALSE,
                                 na.rm=TRUE,
                                 out.rm=NULL,
                                 breakline=FALSE, 
                                 models="default",
                                 cubic=FALSE,
                                 verbose=TRUE,
                                 estimator="MLR",
                                 se="robust",
                                 missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.cognition.relationshipsat)
getPar(RSA.cognition.relationshipsat)
# plot RSA graph 
plot(RSA.cognition.relationshipsat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-7.5,7.5),   
     ylim    = c(-7.5,7.5),
     zlim    = c(1,5),
     xlab = "Cognitive Ability", 
     ylab = "Self View", 
     zlab = "Relationship Satisfaction"
)

#### ROBUSTNESS CHECK WITH COMMON GRAND MEAN CENTERING
#### common grand mean centering
commonmean <- (mean(selfinsight.cog.relationshipsat$raventotal)+mean(selfinsight.cog.relationshipsat$ravenselfrating))/2; commonmean
selfinsight.cog.relationshipsat$cc.raventotal <- (selfinsight.cog.relationshipsat$raventotal)-commonmean
selfinsight.cog.relationshipsat$cc.ravenselfrating <- (selfinsight.cog.relationshipsat$ravenselfrating)-commonmean
with(selfinsight.cog.relationshipsat,describe(cbind(cc.raventotal,raventotal,cc.ravenselfrating,ravenselfrating,commonmean)))
#### test for self-insight with polynomial regression analyses and RSA
RSA.cognition.relationshipsat.cc<-RSA(relationship.satisfaction~cc.raventotal*cc.ravenselfrating,data=selfinsight.cog.relationshipsat,
                                    center=FALSE,
                                    scale=FALSE,
                                    na.rm=TRUE,
                                    out.rm=NULL,
                                    breakline=FALSE, 
                                    models="default",
                                    cubic=FALSE,
                                    verbose=TRUE,
                                    estimator="MLR",
                                    se="robust",
                                    missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.cognition.relationshipsat.cc)
getPar(RSA.cognition.relationshipsat.cc)
# plot RSA graph 
plot(RSA.cognition.relationshipsat.cc,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1"), 
     project = c("PA1", "LOC"),        
     xlim    = c(-10,10),   
     ylim    = c(-10,10),    
     zlim    = c(1,5),
     xlab = "Cognitive Ability", 
     ylab = "Self View", 
     zlab = "Relationship Satisfaction"
)

#### !!!! ALL ANALYSES BELOW ARE EXPLORATORY !!!! ####
#### !!!! NONE OF THE ANALYSES BELOW ARE DONE TO CONDUCT CONFIRMATORY TESTS OF HYPOTHESES !!!!

#### EXPLORATORY CONDITION-BASED REGRESSION ANALYSIS TO TEST SELF-ENHANCEMENT ####
#### code is adapted from Humberg et al. 

## first test: emotional abilities and life satisfaction
## Justify moving from a polynomial model to linear model
summary(full_model.1.linear <- lm(life.satisfaction~c.ERAtotal+c.ertselfrating, data=selfinsight.emo.lifesat, na.action=na.exclude))
summary(full_model.1.poly <- lm(life.satisfaction~c.ERAtotal+c.ertselfrating+c.ERAtotal2+c.ERAtotal_c.ertselfrating+c.ertselfrating2, data=selfinsight.emo.lifesat, na.action=na.exclude))
anova(full_model.1.linear, full_model.1.poly)
## Specify the regression model, where the outcome is predicted by R and S
model.emo.life <- lm(life.satisfaction ~ ERAtotal + ertselfrating, data=selfinsight.emo.lifesat, na.action=na.omit); summary(model.emo.life)
model <- 'life.satisfaction ~ 1 + c1*ERAtotal + c2*ertselfrating
a1 := c1+c2
a3 := c2-c1
c1_times_2 := 2*c1
c2_times_2 := 2*c2
minus_c1_times_2 := -2*c1
minus_c2_times_2 := -2*c2'
## Compute the regression with robust SE
regression <- sem(model, se="robust.huber.white", data=selfinsight.emo.lifesat)
## Table of estimated coefficients
estimates <- parameterEstimates(regression)
## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
## its standard error, respective p-value, and the sign of (c2 - c1)
## with the help of the lemma in Supplemental Material G:
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "c2_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "c1_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
}
print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)

## second test: emotional abilities and career satisfaction
## Justify moving from a polynomial model to linear model
summary(full_model.2.linear <- lm(career.satisfaction~c.ERAtotal+c.ertselfrating, data=selfinsight.emo.careersat, na.action=na.exclude))
summary(full_model.2.poly <- lm(career.satisfaction~c.ERAtotal+c.ertselfrating+c.ERAtotal2+c.ERAtotal_c.ertselfrating+c.ertselfrating2, data=selfinsight.emo.careersat, na.action=na.exclude))
anova(full_model.2.linear, full_model.2.poly)
## Specify the regression model, where the outcome is predicted by R and S
model.emo.career <- lm(career.satisfaction ~ ERAtotal + ertselfrating, data=selfinsight.emo.careersat, na.action=na.omit); summary(model.emo.career)
model <- 'career.satisfaction ~ 1 + c1*ERAtotal + c2*ertselfrating
a1 := c1+c2
a3 := c2-c1
c1_times_2 := 2*c1
c2_times_2 := 2*c2
minus_c1_times_2 := -2*c1
minus_c2_times_2 := -2*c2'
## Compute the regression
regression <- sem(model, se="robust.huber.white", data=selfinsight.emo.careersat)
## Table of estimated coefficients
estimates <- parameterEstimates(regression)
## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
## its standard error, respective p-value, and the sign of (c2 - c1)
## with the help of the lemma in Supplemental Material G:
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "c2_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "c1_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
}
print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)

## third test: emotional abilities and relationship satisfaction
#### Justify moving from a polynomial model to linear model
summary(full_model.3.linear <- lm(relationship.satisfaction~c.ERAtotal+c.ertselfrating, data=selfinsight.emo.relationshipsat, na.action=na.exclude))
summary(full_model.3.poly <- lm(relationship.satisfaction~c.ERAtotal+c.ertselfrating+c.ERAtotal2+c.ERAtotal_c.ertselfrating+c.ertselfrating2, data=selfinsight.emo.relationshipsat, na.action=na.exclude))
anova(full_model.3.linear, full_model.3.poly)
# polynomial model is significantly better; not justified to move to linear model
## exploratory CRA tests not justified; not justified to move from full polynomial model to linear model
# ## Specify the regression model, where the outcome is predicted by R and S
# model.emo.relationship2 <- lm(relationship.satisfaction ~ c.ERAtotal + c.ertselfrating, data=selfinsight.emo.relationshipsat, na.action=na.omit); summary(model.emo.relationship)
# model <- 'relationship.satisfaction ~ 1 + c1*ERAtotal + c2*ertselfrating 
# a1 := c1+c2
# a3 := c2-c1
# c1_times_2 := 2*c1
# c2_times_2 := 2*c2
# minus_c1_times_2 := -2*c1
# minus_c2_times_2 := -2*c2'
# ## Compute the regression
# regression <- sem(model, se="robust.huber.white", data=selfinsight.emo.relationshipsat)
# ## Table of estimated coefficients
# estimates <- parameterEstimates(regression)
# ## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
# ## its standard error, respective p-value, and the sign of (c2 - c1)
# ## with the help of the lemma in Supplemental Material G:
# if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
#   abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
#   se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
#   a3_is <- "positive"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
# }
# if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
#   abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
#   se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
#   a3_is <- "negative"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
# }
# if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
#   abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
#   se <- round(subset(estimates, label == "c2_times_2")["se"],2)
#   a3_is <- "positive"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
# }
# if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
#   abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
#   se <- round(subset(estimates, label == "c1_times_2")["se"],2)
#   a3_is <- "negative"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
# }
# print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)

## fourth test: cognitive abilities and life satisfaction
## Justify moving from a polynomial model to linear model
summary(full_model.4.linear <- lm(life.satisfaction~c.raventotal+c.ravenselfrating, data=selfinsight.cog.lifesat, na.action=na.exclude))
summary(full_model.4.poly <- lm(life.satisfaction~c.raventotal+c.ravenselfrating+c.raventotal2+c.raventotal_c.ravenselfrating+c.ravenselfrating2, data=selfinsight.cog.lifesat, na.action=na.exclude))
anova(full_model.4.linear, full_model.4.poly)
## The overall model is not significant; CRA analysis not appropriate
# ## Specify the regression model, where the outcome is predicted by R and S
# model.cog.life <- lm(life.satisfaction ~ raventotal + ravenselfrating, data=selfinsight.cog.lifesat, na.action=na.omit); summary(model.cog.life)
# model <- 'life.satisfaction ~ 1 + c1*raventotal + c2*ravenselfrating
# a1 := c1+c2
# a3 := c2-c1
# c1_times_2 := 2*c1
# c2_times_2 := 2*c2
# minus_c1_times_2 := -2*c1
# minus_c2_times_2 := -2*c2'
# ## Compute the regression
# regression <- sem(model, se="robust.huber.white", data=selfinsight.cog.lifesat)
# ## Table of estimated coefficients
# estimates <- parameterEstimates(regression)
# ## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
# ## its standard error, respective p-value, and the sign of (c2 - c1)
# ## with the help of the lemma in Supplemental Material G:
# if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
#   abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
#   se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
#   a3_is <- "positive"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
# }
# if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
#   abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
#   se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
#   a3_is <- "negative"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
# }
# if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
#   abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
#   se <- round(subset(estimates, label == "c2_times_2")["se"],2)
#   a3_is <- "positive"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
# }
# if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
#   abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
#   se <- round(subset(estimates, label == "c1_times_2")["se"],2)
#   a3_is <- "negative"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
# }
# print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
# ## Interpretation:
# ## abs is NOT significantly bigger than zero
# ## therefore, there is NO effect of self enhancement (SE)

## fifth test: cognitive abilities and career satisfaction
##Justify moving from a polynomial model to linear model
summary(full_model.5.linear <- lm(career.satisfaction~c.raventotal+c.ravenselfrating, data=selfinsight.cog.careersat, na.action=na.exclude))
summary(full_model.5.poly <- lm(career.satisfaction~c.raventotal+c.ravenselfrating+c.raventotal2+c.raventotal_c.ravenselfrating+c.ravenselfrating2, data=selfinsight.cog.careersat, na.action=na.exclude))
anova(full_model.5.linear, full_model.5.poly)
## Specify the regression model, where the outcome is predicted by R and S
model.cog.career <- lm(career.satisfaction ~ raventotal + ravenselfrating, data=selfinsight.cog.careersat, na.action=na.omit); summary(model.cog.career)
model <- 'career.satisfaction ~ 1 + c1*raventotal + c2*ravenselfrating
a1 := c1+c2
a3 := c2-c1
c1_times_2 := 2*c1
c2_times_2 := 2*c2
minus_c1_times_2 := -2*c1
minus_c2_times_2 := -2*c2'
## Compute the regression
regression <- sem(model, se="robust.huber.white", data=selfinsight.cog.careersat)
## Table of estimated coefficients
estimates <- parameterEstimates(regression)
## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
## its standard error, respective p-value, and the sign of (c2 - c1)
## with the help of the lemma in Supplemental Material G:
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "c2_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "c1_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
}
print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)

## sixth test: cognitive abilities and relationship satisfaction
## Justify moving from a polynomial model to linear model
summary(full_model.6.linear <- lm(relationship.satisfaction~c.raventotal+c.ravenselfrating, data=selfinsight.cog.relationshipsat, na.action=na.exclude))
summary(full_model.6.poly <- lm(relationship.satisfaction~c.raventotal+c.ravenselfrating+c.raventotal2+c.raventotal_c.ravenselfrating+c.ravenselfrating2, data=selfinsight.cog.relationshipsat, na.action=na.exclude))
anova(full_model.6.linear, full_model.6.poly)
## Specify the regression model, where the outcome is predicted by R and S
model.cog.relationship <- lm(relationship.satisfaction ~ raventotal+  ravenselfrating, data=selfinsight.cog.relationshipsat, na.action=na.omit); summary(model.cog.relationship)
model <- 'relationship.satisfaction ~ 1 + c1*raventotal + c2*ravenselfrating 
a1 := c1+c2
a3 := c2-c1
c1_times_2 := 2*c1
c2_times_2 := 2*c2
minus_c1_times_2 := -2*c1
minus_c2_times_2 := -2*c2'
## Compute the regression
regression <- sem(model, se="robust.huber.white", data=selfinsight.cog.relationshipsat)
## Table of estimated coefficients
estimates <- parameterEstimates(regression)
## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
## its standard error, respective p-value, and the sign of (c2 - c1)
## with the help of the lemma in Supplemental Material G:
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "c2_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "c1_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
}
print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)

#### EXPLORATORY ANALYSES WITH STANDARDIZED PREDICTORS ####

#### emotional abilities and life satisfaction
selfinsight.emo.lifesat$z.ERAtotal <- (selfinsight.emo.lifesat$ERAtotal-mean(selfinsight.emo.lifesat$ERAtotal,na.rm=TRUE))/sd(selfinsight.emo.lifesat$ERAtotal,na.rm=TRUE)
selfinsight.emo.lifesat$z.ertselfrating <- (selfinsight.emo.lifesat$ertselfrating-mean(selfinsight.emo.lifesat$ertselfrating,na.rm=TRUE))/sd(selfinsight.emo.lifesat$ertselfrating,na.rm=TRUE)
with(selfinsight.emo.lifesat, describe(cbind(z.ERAtotal,z.ertselfrating,ERAtotal,ertselfrating)))
RSA.STD.emotion.lifesat<-RSA(life.satisfaction~z.ERAtotal*z.ertselfrating,data=selfinsight.emo.lifesat,
                         center=FALSE,
                         scale=FALSE,
                         na.rm=TRUE,
                         out.rm=NULL,
                         breakline=FALSE, 
                         models="default",
                         cubic=FALSE,
                         verbose=TRUE,
                         estimator="MLR",
                         se="robust",
                         missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.STD.emotion.lifesat)
getPar(RSA.STD.emotion.lifesat)
# plot RSA graph 
plot(RSA.STD.emotion.lifesat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1", "PA2"), 
     project = c("PA1", "PA2", "LOC"),        
     xlim    = c(-5,5),   
     ylim    = c(-5,5),    
     xlab = "STD Emotional Ability", 
     ylab = "STD Self View", 
     zlab = "Life Satisfaction"
)
## Specify the regression model, where the outcome is predicted by R and S
model.emo.life <- lm(life.satisfaction ~  z.ERAtotal + z.ertselfrating, data=selfinsight.emo.lifesat, na.action=na.omit); summary(model.emo.life)
model <- 'life.satisfaction ~ 1 + c1*z.ERAtotal + c2*z.ertselfrating
a1 := c1+c2
a3 := c2-c1
c1_times_2 := 2*c1
c2_times_2 := 2*c2
minus_c1_times_2 := -2*c1
minus_c2_times_2 := -2*c2'
## Compute the regression
regression <- sem(model, se="robust.huber.white", data=selfinsight.emo.lifesat)
## Table of estimated coefficients
estimates <- parameterEstimates(regression)
## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
## its standard error, respective p-value, and the sign of (c2 - c1)
## with the help of the lemma in Supplemental Material G:
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "c2_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "c1_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
}
print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)


## power analysis
model.emo.life <- lm(life.satisfaction ~ z.ERAtotal, data=selfinsight.emo.lifesat, na.action=na.omit); summary(model.emo.life)
model.emo.life <- lm(life.satisfaction ~ z.ertselfrating, data=selfinsight.emo.lifesat, na.action=na.omit); summary(model.emo.life)
model.emo.life <- lm(life.satisfaction ~ z.ertselfrating + z.ERAtotal, data=selfinsight.emo.lifesat, na.action=na.omit); summary(model.emo.life)


#### emotional abilities and career satisfaction
selfinsight.emo.careersat$z.ERAtotal <- (selfinsight.emo.careersat$ERAtotal-mean(selfinsight.emo.careersat$ERAtotal,na.rm=TRUE))/sd(selfinsight.emo.careersat$ERAtotal,na.rm=TRUE)
selfinsight.emo.careersat$z.ertselfrating <- (selfinsight.emo.careersat$ertselfrating-mean(selfinsight.emo.careersat$ertselfrating,na.rm=TRUE))/sd(selfinsight.emo.careersat$ertselfrating,na.rm=TRUE)
with(selfinsight.emo.careersat, describe(cbind(z.ERAtotal,z.ertselfrating,ERAtotal,ertselfrating)))
RSA.STD.emotion.careersat<-RSA(career.satisfaction~z.ERAtotal*z.ertselfrating,data=selfinsight.emo.careersat,
                             center=FALSE,
                             scale=FALSE,
                             na.rm=TRUE,
                             out.rm=NULL,
                             breakline=FALSE, 
                             models="default",
                             cubic=FALSE,
                             verbose=TRUE,
                             estimator="MLR",
                             se="robust",
                             missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.STD.emotion.careersat)
getPar(RSA.STD.emotion.careersat)
# plot RSA graph 
plot(RSA.STD.emotion.careersat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1", "PA2"), 
     project = c("PA1", "PA2", "LOC"),        
     xlim    = c(-5,5),   
     ylim    = c(-5,5),    
     xlab = "STD Emotional Ability", 
     ylab = "STD Self View", 
     zlab = "career Satisfaction"
)
## Specify the regression model, where the outcome is predicted by R and S
model.emo.career <- lm(career.satisfaction ~  z.ERAtotal + z.ertselfrating, data=selfinsight.emo.careersat, na.action=na.omit); summary(model.emo.career)
model <- 'career.satisfaction ~ 1 + c1*z.ERAtotal + c2*z.ertselfrating
a1 := c1+c2
a3 := c2-c1
c1_times_2 := 2*c1
c2_times_2 := 2*c2
minus_c1_times_2 := -2*c1
minus_c2_times_2 := -2*c2'
## Compute the regression
regression <- sem(model, se="robust.huber.white", data=selfinsight.emo.careersat)
## Table of estimated coefficients
estimates <- parameterEstimates(regression)
## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
## its standard error, respective p-value, and the sign of (c2 - c1)
## with the help of the lemma in Supplemental Material G:
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "c2_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "c1_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
}
print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)

#### emotional abilities and relationship satisfaction
selfinsight.emo.relationshipsat$z.ERAtotal <- (selfinsight.emo.relationshipsat$ERAtotal-mean(selfinsight.emo.relationshipsat$ERAtotal,na.rm=TRUE))/sd(selfinsight.emo.relationshipsat$ERAtotal,na.rm=TRUE)
selfinsight.emo.relationshipsat$z.ertselfrating <- (selfinsight.emo.relationshipsat$ertselfrating-mean(selfinsight.emo.relationshipsat$ertselfrating,na.rm=TRUE))/sd(selfinsight.emo.relationshipsat$ertselfrating,na.rm=TRUE)
with(selfinsight.emo.relationshipsat, describe(cbind(z.ERAtotal,z.ertselfrating,ERAtotal,ertselfrating)))
RSA.STD.emotion.relationshipsat<-RSA(relationship.satisfaction~z.ERAtotal*z.ertselfrating,data=selfinsight.emo.relationshipsat,
                               center=FALSE,
                               scale=FALSE,
                               na.rm=TRUE,
                               out.rm=NULL,
                               breakline=FALSE, 
                               models="default",
                               cubic=FALSE,
                               verbose=TRUE,
                               estimator="MLR",
                               se="robust",
                               missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.STD.emotion.relationshipsat)
getPar(RSA.STD.emotion.relationshipsat)
# plot RSA graph 
plot(RSA.STD.emotion.relationshipsat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1", "PA2"), 
     project = c("PA1", "PA2", "LOC"),        
     xlim    = c(-5,5),   
     ylim    = c(-5,5),    
     xlab = "STD Emotional Ability", 
     ylab = "STD Self View", 
     zlab = "relationship Satisfaction"
)
## Specify the regression model, where the outcome is predicted by R and S
model.emo.relationship <- lm(relationship.satisfaction ~  z.ERAtotal + z.ertselfrating, data=selfinsight.emo.relationshipsat, na.action=na.omit); summary(model.emo.relationship)
model <- 'relationship.satisfaction ~ 1 + c1*z.ERAtotal + c2*z.ertselfrating
a1 := c1+c2
a3 := c2-c1
c1_times_2 := 2*c1
c2_times_2 := 2*c2
minus_c1_times_2 := -2*c1
minus_c2_times_2 := -2*c2'
## Compute the regression
regression <- sem(model, se="robust.huber.white", data=selfinsight.emo.relationshipsat)
## Table of estimated coefficients
estimates <- parameterEstimates(regression)
## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
## its standard error, respective p-value, and the sign of (c2 - c1)
## with the help of the lemma in Supplemental Material G:
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "c2_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "c1_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
}
print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)

#### cognitive abilities and life satisfaction
selfinsight.cog.lifesat$z.raventotal <- (selfinsight.cog.lifesat$raventotal-mean(selfinsight.cog.lifesat$raventotal,na.rm=TRUE))/sd(selfinsight.cog.lifesat$raventotal,na.rm=TRUE)
selfinsight.cog.lifesat$z.ravenselfrating <- (selfinsight.cog.lifesat$ravenselfrating-mean(selfinsight.cog.lifesat$ravenselfrating,na.rm=TRUE))/sd(selfinsight.cog.lifesat$ravenselfrating,na.rm=TRUE)
with(selfinsight.cog.lifesat, describe(cbind(z.raventotal,z.ravenselfrating,raventotal,ravenselfrating)))
RSA.STD.emotion.lifesat<-RSA(life.satisfaction~z.raventotal*z.ravenselfrating,data=selfinsight.cog.lifesat,
                               center=FALSE,
                               scale=FALSE,
                               na.rm=TRUE,
                               out.rm=NULL,
                               breakline=FALSE, 
                               models="default",
                               cubic=FALSE,
                               verbose=TRUE,
                               estimator="MLR",
                               se="robust",
                               missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.STD.emotion.lifesat)
getPar(RSA.STD.emotion.lifesat)
# plot RSA graph 
plot(RSA.STD.emotion.lifesat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1", "PA2"), 
     project = c("PA1", "PA2", "LOC"),        
     xlim    = c(-3,3),   
     ylim    = c(-3,3),    
     xlab = "STD cognitive Ability", 
     ylab = "STD Self View", 
     zlab = "Life Satisfaction"
)
## Specify the regression model, where the outcome is predicted by R and S
## Linear model is not significant; CRA analyses not conducted
# model.cog.life <- lm(life.satisfaction ~  z.raventotal + z.ravenselfrating, data=selfinsight.cog.lifesat, na.action=na.omit); summary(model.cog.life)
# model <- 'life.satisfaction ~ 1 + c1*z.raventotal + c2*z.ravenselfrating
# a1 := c1+c2
# a3 := c2-c1
# c1_times_2 := 2*c1
# c2_times_2 := 2*c2
# minus_c1_times_2 := -2*c1
# minus_c2_times_2 := -2*c2'
# ## Compute the regression
# regression <- sem(model, se="robust.huber.white", data=selfinsight.cog.lifesat)
# ## Table of estimated coefficients
# estimates <- parameterEstimates(regression)
# ## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
# ## its standard error, respective p-value, and the sign of (c2 - c1)
# ## with the help of the lemma in Supplemental Material G:
# if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
#   abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
#   se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
#   a3_is <- "positive"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
# }
# if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
#   abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
#   se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
#   a3_is <- "negative"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
# }
# if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
#   abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
#   se <- round(subset(estimates, label == "c2_times_2")["se"],2)
#   a3_is <- "positive"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
# }
# if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
#   abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
#   se <- round(subset(estimates, label == "c1_times_2")["se"],2)
#   a3_is <- "negative"
#   # compute one-tailed p-value of abs for the hypothesis abs > 0, 
#   # depending on the tail in which abs is positioned
#   if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
#   if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
# }
# print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
# ## Interpretation:
# ## abs is NOT significantly bigger than zero
# ## therefore, there is NO effect of self enhancement (SE)

#### cognitive abilities and career satisfaction
selfinsight.cog.careersat$z.raventotal <- (selfinsight.cog.careersat$raventotal-mean(selfinsight.cog.careersat$raventotal,na.rm=TRUE))/sd(selfinsight.cog.careersat$raventotal,na.rm=TRUE)
selfinsight.cog.careersat$z.ravenselfrating <- (selfinsight.cog.careersat$ravenselfrating-mean(selfinsight.cog.careersat$ravenselfrating,na.rm=TRUE))/sd(selfinsight.cog.careersat$ravenselfrating,na.rm=TRUE)
with(selfinsight.cog.careersat, describe(cbind(z.raventotal,z.ravenselfrating,raventotal,ravenselfrating)))
RSA.STD.emotion.careersat<-RSA(career.satisfaction~z.raventotal*z.ravenselfrating,data=selfinsight.cog.careersat,
                               center=FALSE,
                               scale=FALSE,
                               na.rm=TRUE,
                               out.rm=NULL,
                               breakline=FALSE, 
                               models="default",
                               cubic=FALSE,
                               verbose=TRUE,
                               estimator="MLR",
                               se="robust",
                               missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.STD.emotion.careersat)
getPar(RSA.STD.emotion.careersat)
# plot RSA graph 
plot(RSA.STD.emotion.careersat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1", "PA2"), 
     project = c("PA1", "PA2", "LOC"),        
     xlim    = c(-3,3),   
     ylim    = c(-3,3),    
     xlab = "STD cognitive Ability", 
     ylab = "STD Self View", 
     zlab = "Career Satisfaction"
)
## Specify the regression model, where the outcome is predicted by R and S
model.cog.career <- lm(career.satisfaction ~  z.raventotal + z.ravenselfrating, data=selfinsight.cog.careersat, na.action=na.omit); summary(model.cog.career)
model <- 'career.satisfaction ~ 1 + c1*z.raventotal + c2*z.ravenselfrating
a1 := c1+c2
a3 := c2-c1
c1_times_2 := 2*c1
c2_times_2 := 2*c2
minus_c1_times_2 := -2*c1
minus_c2_times_2 := -2*c2'
## Compute the regression
regression <- sem(model, se="robust.huber.white", data=selfinsight.cog.careersat)
## Table of estimated coefficients
estimates <- parameterEstimates(regression)
## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
## its standard error, respective p-value, and the sign of (c2 - c1)
## with the help of the lemma in Supplemental Material G:
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "c2_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "c1_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
}
print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)

#### cognitive abilities and relationship satisfaction
selfinsight.cog.relationshipsat$z.raventotal <- (selfinsight.cog.relationshipsat$raventotal-mean(selfinsight.cog.relationshipsat$raventotal,na.rm=TRUE))/sd(selfinsight.cog.relationshipsat$raventotal,na.rm=TRUE)
selfinsight.cog.relationshipsat$z.ravenselfrating <- (selfinsight.cog.relationshipsat$ravenselfrating-mean(selfinsight.cog.relationshipsat$ravenselfrating,na.rm=TRUE))/sd(selfinsight.cog.relationshipsat$ravenselfrating,na.rm=TRUE)
with(selfinsight.cog.relationshipsat, describe(cbind(z.raventotal,z.ravenselfrating,raventotal,ravenselfrating)))
RSA.STD.emotion.relationshipsat<-RSA(relationship.satisfaction~z.raventotal*z.ravenselfrating,data=selfinsight.cog.relationshipsat,
                                     center=FALSE,
                                     scale=FALSE,
                                     na.rm=TRUE,
                                     out.rm=NULL,
                                     breakline=FALSE, 
                                     models="default",
                                     cubic=FALSE,
                                     verbose=TRUE,
                                     estimator="MLR",
                                     se="robust",
                                     missing='listwise'
)
# see summary and parameters of RSA
summary(RSA.STD.emotion.relationshipsat)
getPar(RSA.STD.emotion.relationshipsat)
# plot RSA graph 
plot(RSA.STD.emotion.relationshipsat,      
     main = NULL,      
     param   = FALSE,     
     coefs   = FALSE,      
     axes    = c("LOC", "LOIC", "PA1", "PA2"), 
     project = c("PA1", "PA2", "LOC"),        
     xlim    = c(-5,5),   
     ylim    = c(-5,5),    
     xlab = "STD cognitive Ability", 
     ylab = "STD Self View", 
     zlab = "Relationship Satisfaction"
)
## Specify the regression model, where the outcome is predicted by R and S
model.cog.relationship <- lm(relationship.satisfaction ~  z.raventotal + z.ravenselfrating, data=selfinsight.cog.relationshipsat, na.action=na.omit); summary(model.cog.relationship)
model <- 'relationship.satisfaction ~ 1 + c1*z.raventotal + c2*z.ravenselfrating
a1 := c1+c2
a3 := c2-c1
c1_times_2 := 2*c1
c2_times_2 := 2*c2
minus_c1_times_2 := -2*c1
minus_c2_times_2 := -2*c2'
## Compute the regression
regression <- sem(model, se="robust.huber.white", data=selfinsight.cog.relationshipsat)
## Table of estimated coefficients
estimates <- parameterEstimates(regression)
## Determination of the parameter abs := |c2 - c1| - |c1 + c2|,
## its standard error, respective p-value, and the sign of (c2 - c1)
## with the help of the lemma in Supplemental Material G:
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c1_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c1_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]>0){
  abs <- round(subset(estimates, label == "minus_c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "minus_c2_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "minus_c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "minus_c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]>=0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c2_times_2")["est"],2)
  se <- round(subset(estimates, label == "c2_times_2")["se"],2)
  a3_is <- "positive"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c2_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c2_times_2")["pvalue"]/2), 5)}
}
if (subset(estimates, label == "a3")["est"]<0 & subset(estimates, label == "a1")["est"]<0){
  abs <- round(subset(estimates, label == "c1_times_2")["est"],2)
  se <- round(subset(estimates, label == "c1_times_2")["se"],2)
  a3_is <- "negative"
  # compute one-tailed p-value of abs for the hypothesis abs > 0, 
  # depending on the tail in which abs is positioned
  if (abs >= 0){pvalue <- round(subset(estimates, label == "c1_times_2")["pvalue"]/2, 5)}
  if (abs < 0){pvalue <- round(1 - (subset(estimates, label == "c1_times_2")["pvalue"]/2), 5)}
}
print(c("abs"=abs,se,pvalue,"a3 is"=a3_is))
## Interpretation:
## abs is significantly bigger than zero
## therefore, there is an effect of self enhancement (SE) (abs > 0) 
## this effect has the same sign as a3 := (c2 - c1). the sign is is positive (c2 - c1 > 0)

#### EXPLORE GENDER DIFFERENCE ####

#### correlate gender (men and women) with abilities test scores with self-views of abilities
selfinsight.menwomen <- subset(selfinsight, selfinsight$gendercoded==1|selfinsight$gendercoded==2)
with(selfinsight.menwomen, corr.test(cbind(gendercoded,ERAtotal,ertselfrating,raventotal,ravenselfrating)), use="pairwise", adjust="none")

